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ABSTRACT 


A  great  deal  of  controversy  exists  concerning  the  possible  global  wanning  trend 
which  may  occur  as  a  result  of  a  documented  increase  in  atmospheric  "greenhouse"  gasses. 
The  1991  Heard  Island  Feasibility  Experiment  tested  the  feasibility  of  using  transmissions 
of  acoustic  energy  through  the  major  ocean  basins  of  the  world  to  monitor  spatially 
averaged  global  temperature  trends.  This  thesis  documents  the  Naval  Postgraduate 
School's  reception  of  the  phase  encoded  signal  transmitted  from  the  Southern  Indian 
Ocean,  development  of  real-time  signal  processing  software,  and  preliminary  data  analysis. 

Data,  received  from  a  32-channel  vertical  array  suspended  in  the  deep  sound  channel 
off  the  coast  of  Monterey,  CA,  was  processed  using  real-time  capable  software.  Data 
processing  to  reduce  noise,  determine  SNR,  and  remove  the  m- sequence  coding  was  found 
to  be  quite  sensitive  to  Doppler  frequency  shifts.  Although  the  SNR  of  the  raw  data  was 
only  about  -27.5  dB  for  individual  hydrophones,  the  transmitted  signal  was  detected  in 
both  the  frequency  and  time  domains.  However,  the  maximum  processed  signal  peak  in 
the  time  domain  had  an  SNR  of  only  +9  dB  which  is  insufficient  for  use  in  a  long  term 
global  temperature  monitoring  project,  as  it  provides  inadequate  arrival  time  resolution. 
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I.  INTRODUCTION 


A.  GLOBAL  ACOUSTIC  TRANSMISSIONS 

The  widely  discussed  theory  of  greenhouse  wanning  suggests  that  increasing  levels 
of  certain  atmospheric  gasses,  such  as  carbon  dioxide  and  methane,  will  ultimately  lead 
to  climatological  wanning  on  a  global  scale.  Well  documented  increases  in  atmospheric 
carbon  dioxide  levels  from  air  (Keeling  and  others,  1989)  and  ice  core  samples  (Neftel 
and  others,  1985)  are  shown  in  Figure  1.1. 
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Studies  of  global  temperature  trends,  including  extensive  numerical  modeling, 
indicate  a  possible  correlation  between  the  increase  of  these  gasses  and  planetary 
wanning.  In  one  study  (Hansen  and  Lebedeff,  1987)  the  historical  atmospheric  surface 
temperature  records  over  the  last  century  were  compiled  with  the  averaged  results 
displayed  in  Figure  1.2.  In  addition,  the  Carbon  Dioxide  Assessment  Committee  of  the 
National  Research  Council  has  surveyed  the  results  of  numerous  global  circulation  models 
to  project  a  continued  warming  trend  ranging  from  1  to  5  degrees  Celsius  within  the  next 
50  years  (NRC,  1983).  This  magnitude  of  temperature  change  could  have  dramatic 
consequences  for  the  world  ecosystem  with  significant  environmental,  economic,  and 


1880  1900  1920  1940  1960  1980 

Figure  1.2:  Global  Temperature  Trend  (Hansen  and  Lebedeff,  1987). 
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political  ramifications.  Understandably,  greenhouse  warming  is  the  subject  of 
considerable  scientific  and  political  debate,  as  numerous  challenging  questions  remain 
unanswered  concerning  the  magnitude,  localization,  and  rate  of  wanning.  However,  since 
addressing  the  problem  of  global  warming  could  involve  huge  expenditures  of  human 
resources,  prudent  planning  and  allocation  of  limited  assets  requires  a  more  precise 
understanding  of  the  coupling  between  the  "greenhouse"  gas  concentration  and  planetary 
temperatures. 

A  reconstruction  of  global  air  temperature  records  over  the  past  108  years,  shown 
in  Figure  1.2,  indicates  an  overall  temperature  increase  of  0.5  °C  and  a  1990  warming  rate 
of  0.02  °C  per  year.  However,  in  order  to  statistically  detect  a  trend  in  observed 
temperature  data  at  a  95%  confidence  level,  a  time  series  long  enough  to  produce  a 
change  in  the  actual  mean  of  four  times  the  standard  deviation,  <j,  is  required  (Munk  and 
Forbes,  1989).  The  Keeling  C02  data  shown  in  Figure  1.1  clearly  shows  this  long  term 
trend  with  over  40  a  over  33  years.  Due  to  the  inherent  variability  of  atmospheric 
temperatures,  the  4  a  statistical  criteria  would  require  approximately  100  years  of  air 
temperature  measurements  to  obtain  a  95  percent  certainty  that  the  earth’s  climate  is 
indeed  changing. 

Munk  and  Forbes  have  suggested  obtaining  a  time  series  of  integrated  ocean 
temperatures  which,  because  of  their  smaller  variance  and  direct  coupling  with  the 
atmosphere,  could  produce  the  same  degree  of  certainty  in  only  10  years.  They  have 
proposed  measuring  the  increase  in  sound  speed  through  all  major  ocean  basins  using 
global  acoustic  transmissions  of  low  frequency  sound.  The  highly  efficient  waveguide 
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which  results  from  the  oceanic  sound  speed  minimum  near  a  depth  of  1  km  at  lower  and 
mid-latitudes  gives  reasonable  expectation  of  global  sound  propagation  from  a  strong 
acoustic  source.  Since  temperature  is  the  dominant  influence  on  the  speed  of  sound  in 
seawater  (increasing  by  5  m/s/°C),  the  resulting  travel  time  change  would  provide  a 
measure  of  the  mean  temperature  change  along  the  acoustic  path. 

Predicted  oceanic  wanning  of  0.02  °C/year  at  the  surface  and  0.004  °C  at  the  sound 
channel  axis  would  decrease  the  travel  time  over  a  15,000  km  range  by  0.25  seconds  per 
year  (Munk  and  Forbes,  1989).  This  reduction  of  travel  time  w  ild  result  in  a  4  a 
change  (2.5  °C)  in  about  10  years.  Figure  1,3  shows  the  expected  decrease  in  travel  time 
from  the  Princeton  model  of  greenhouse  warming  (straight  line)  and  a  multi-level 
primitive-equation  mode  1  output  of  mesoscale  variability  (Manabe  and  Stouffer,  1986  and 
Semtner  and  Chervin,  1990).  The  planetary  scale  of  the  transmissions  would  smooth  out 
any  localized  temperature  variations  associated  with  mesoscale  eddies  to  provide 


Figure  1.3:  Signal  Travel  Time  vs  Calendar  Years  (Munk  and  Forbes,  1989). 
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indications  of  gyre  and  basin  scale  wanning.  The  average  travel  time  of  the  acoustic 
signal  would  thus  provide  an  indication  of  the  integrated  temperature  for  each  travel  path. 


B.  THE  1991  HEARD  ISLAND  FEASIBILITY  EXPERIMENT 


1.  Objectives 

In  January  of  1991  scientists  representing  9  nations  and  13  universities  and 
research  institutions  embarked  on  a  worldwide  experiment  to  test  the  feasibility  of  using 
acoustic  transmissions  to  monitor  the  average  planetary  temperature.  An  array  of 
powerful  acoustic  transducers  was  suspended  in  the  sound  channel  approximately  40  miles 
east  of  Heard  Island  in  the  Southern  Indian  Ocean.  The  total  output  power  of  the  array 
was  approximately  213  dB  referenced  to  1  jiPa  at  1  meter.  Hydrophone  receivers  as 
shown  in  Figure  1 .4  were  deployed  at  the  following  sites: 

INSTITUTE.  COUNTRY 
NOAA,  USA 
U.Auk.,  N. Zealand 
CSIRO,  Australia 
CSIRO,  Australia 
CSIRO,  Australia 
INSU-TAAF,  France 
NIO,  India 
South  Africa 
Japan 

I OS,  Canada 
MBARI/NPS/WHOI/ 
MIT,  USA 
U.  Miami,  USA 
SAIC,  USA 
U.  Mich.,  USA 
U.  Mich.,  USA 


TYPE 


LOCATION 

Ascension  Island 
Tasman  Sea 
Tasmania 

Christmas  Island 
Mawson  Station 
Kerguelen  Island 
Indian  Ocean 
Capetown 
Fiji 

North  Pacific 
Monterey,  CA 

Florida 

Bermuda 

Atlantic  Ocean 
Pacific  Ocean 


Bottom 

Towed  &  Sonobuoy 

Sonobuoys 

Sonobuoys 

Sonobuoys 

Sonobuoys 

Sonobuoys 

Bottom 

Sonobuoys 

Towed  Array 

Vertical  Array 

Swade  Array 
Vertical  Array 
Navy 
Navy 


The  experiment  addressed  three  major  uncertainties  underlying  the  use  of 
global  acoustic  transmissions: 
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If  r  4D*  a*  r  sem  #r  *tr  *tr  tor  ixr  wo*  iw  mo*  iw  f<o*  «o*  *ocr 

Figure  1.4:  Axially  Refracted  Geodesics  (drawn  every  10°  from  the  source)  (Munk 


1990). 

1 .  The  transmitted  signal  strength  required  to  produce  adequate  signal-to-noise  ratio 
at  the  receiving  stations,  including  the  achievable  accuracy  of  arrival  time  measurement, 

2.  The  stability  of  the  multipath  arrival  pattern,  including  the  maximum  coherent 
integration  times,  and 

3.  The  optimum  receiver  locations  for  a  possible  permanent  follow-on  installation. 

2.  Importance 

Successful  conclusion  of  the  Heard  Island  Feasibility  Experiment  entails 
demonstration  of  reliable  propagation  of  low  frequency  acoustic  energy  along  global-scale 
transmission  paths.  Fruitful  results  could  lead  to  the  installation  of  a  permanent  climate 
monitoring  system  which  could  provide  early  warning  of  upper  ocean  warming  of  the 
Southern  Ocean  and  interior  warming  in  other  oceans.  Such  a  capability  could  provide 
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definitive  data  on  the  question  of  global  warming  with  apparent  scientific,  social,  and 
economic  value. 

Since  the  travel  time  of  a  basin-scale  transmission  is  directly  influenced  by 
the  water  characteristics  along  its  path,  the  arrival  structure  would  contain  information  on 
the  meso  and  gyre-scale  processes  through  which  the  signal  transited.  A  permanent 
system  could  also  provide  statistical  information  on  the  variability  of  these  processes 
which  could  be  used  for  calibration  of  operational  eddy  resolving  global  circulation 
models  to  enhance  their  reliability. 

Tactical  employment  of  the  marine  environment  requires  a  broadened 
understanding  of  the  effects  of  meso-scale  effects  on  long-range,  low-frequency  acoustic 
propagation.  Maintaining  the  tactical  edge  in  an  increasingly  transparent  ocean  against 
targets  with  decreasing  acoustic  signatures  demands  exploitation  of  every  available  feature 
to  optimum  advantage.  The  spatial  and  temporal  variability  information  contained  in  the 
arrival  pattern  would  enhance  the  tactician’s  effort  to  optimally  employ  his  vessel.  In 
addition,  the  characteristics  of  a  low-frequency  three-dimensional  wavefield  which  has 
propagated  through  meso-scale  features  will  provide  valuable  information  for  effective 
development  of  future  ASW  systems. 

C.  THESIS  OBJECTIVES 

This  thesis  describes  the  United  States  Naval  Postgraduate  School’s  participation 
in  the  1991  Heard  Island  Feasibility  Experiment  with  an  emphasis  on  signal  processing 
and  preliminary  data  analysis.  The  intent  of  this  thesis  was  threefold: 
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1.  The  first  goal  was  to  deploy  a  vertical  acoustic  array  and  data  collection  hardware 
during  the  at-sea  reception  phase  of  The  Heard  Island  Feasibility  Experiment  to 
digitize  and  store  the  received  signal. 

2.  The  second  objective  was  to  develop  and  implement  signal  processing  software 
modules  in  the  C  programming  language  which  will  interface  with  Laboratory 
Workbench  and  the  Concurrent  6400  real-time  Unix  data  acquisition  computer. 
Laboratory  Workbench  is  a  proprietary  signal  processing  software  product  of  the 
Concurrent  Computer  Corporation. 

3.  The  third  purpose  of  this  thesis  was  to  perform  preliminary  data  analysis  on  the 
data  collected  during  the  at-sea  phase  including  arrival  time  accuracy  estimation, 
determination  of  signal  to  noise  ratio,  and  arrival  structure  analysis. 

D.  THESIS  ORGANIZATION 

The  organization  of  the  remainder  of  this  thesis  is  as  follows: 

Chapter  II  contains  background  information  on  the  1991  Heard  Island  Feasibility 
Experiment,  ocean  acoustic  wave  propagation,  maximal-length  shift-register  sequences, 
Hadamard  transforms,  and  acoustic  Doppler  processing. 

Chapter  HI  specifies  the  data  collection  techniques  and  signal  processing  utilized 
to  analyze  the  received  signal,  including  filtering,  demodulation,  M-sequence  removal, 
data  storage,  and  real  time  processing  considerations. 

Chapter  IV  details  the  results  of  the  data  collection  phase,  the  software  development 
effort,  and  the  results  of  preliminary  data  analysis,  including  resolution,  signal  to  noise 
ratio,  maximum  coherent  integration  time,  multi-path  arrival,  and  depth  profile. 

Chapter  V  summarizes  the  experimental  results,  presents  the  thesis  conclusions  and 
suggests  some  further  efforts  in  analyzing  the  collected  data. 
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II.  BACKGROUND 


A.  THE  1991  HEARD  ISLAND  FEASIBILITY  EXPERIMENT 

1.  Scope  and  Objectives 

The  objectives  of  the  Heard  Island  Feasibility  Experiment  were: 

1.  to  determine  the  required  transmitted  signal  strength  of  a  permanent  installation 
which  would  be  sufficient  to  overcome  propagation  losses  and  provide  an  adequate 
signal-to-noise  ratio  at  the  receiving  stations,  and 

2.  to  determine  whether  the  arrival  pattern  has  features  which  are  sharp  and  stable 
enough  to  determine  the  variable  transmission  time  within  the  10  millisecond 
accuracy  required  to  make  statistically  useful  measurements  of  global  temperature 
trends. 


The  accuracy  of  ocean  basin  average  temperature  estimates  derived  from 
acoustic  propagation  is  dependent  on  the  resolution  of  the  signal’s  arrival  time.  The 
arrival  time  resolution  is  in  turn  related  to  the  bandwidth  of  the  signal  and  the  received 
signal-to-noise  ratio.  Van  Trees  discusses  the  calculation  of  this  type  of  non-linear 
estimate  with  white  Gaussian  noise  (Van  Trees,  1968).  The  result  is  given  as 


°c 


l 

B/SnP. 


(2.1) 


with  a,  the  arrival  time  uncertainty,  B  the  bandwidth  of  the  signal,  and  SNR  the  signal-to- 
noise  ratio  (Spindel,  1985).  Re-arranging  equation  2.1,  the  SNR  required  to  determine 
the  arrival  time  within  o,  is  given  by 


Assuming  a  signal  bandwidth  of  11.4  Hz,  the  desired  arrival  time  accuracy  of  10 
milliseconds  would  require  a  post-processed  receiver  SNR  of  at  least  18.86  dB.  The 
Monterey  effort  of  the  1991  Heard  Island  Experiment  attempted  to  maximize  the  receiver 
SNR  by  optimizing  receiver  locations,  using  phase-encoded  pulse  compression  techniques, 
and  by  using  a  multi-element  vertical  receiving  array. 

2.  Transmission 

It  is  desirable  that  the  permanent  source  for  long  term  global  transmissions 
be  operated  at  less  than  200  dB  rel  pPa@lm  to  allow  reliable  operation  over  a  long 
period  of  time.  For  the  1991  Heard  Island  Feasibility  Experiment,  a  vertical  transmitting 
array  of  ten  Hydroacoustics  HLF4LL  transducers  with  12.5  foot  spacing  was  deployed 
from  the  R/V  Cory  Chouest  between  26  January  and  1  February  1991  at  53°  15’  S,  and 
73°  40’  E,  in  an  area  approximately  40  miles  east  of  Heard  Island  in  the  Southern  Indian 
Ocean.  The  transmitter  depth  was  operationally  limited  to  300  meters  which  mandated 
high  latitude  positioning  where  the  sound  axis  is  shallower.  The  position  near  Heard 
Island  also  allowed  propagation  through  all  major  ocean  basins.  The  array  was  deployed 
such  that  its  center  point  coincided  with  the  deep  sound  channel  axis  near  200  meters. 
The  vessel  was  operated  for  minimal  acceleration,  generating  minimum  thrust  to  maintain 
steerage,  thus  introducing  linear  vice  second  order  Doppler  components  into  the  signal. 
The  operating  radius  of  the  R/V  Cory  Chouest  was  about  15  nautical  miles  which  d;d  not 


affect  the  experiment  since  the  objective  was  to  measure  the  accuracy  of  arrival  time  and 
not  the  time  of  arrival  itself. 

Under  ideal  conditions,  the  total  output  power  of  a  linear  array  of  projectors 

is  given  by 


*te oun  *  Po  +  10  log  (  Nr  ) 


(2.3) 


where  NT  is  the  number  of  active  projectors  and  Pc  is  the  output  power  of  an  individual 
projector.  Nominally,  five  transducers,  alternately  spaced,  were  operated  simultaneously 
and  in  phase  to  produce  a  broadside  beam.  The  operating  transducers  were  changed  as 
required  to  prevent  overheating.  Each  of  the  five  active  transducers  optimally  produced 
206  dB  rel  pPa@lm  to  produce  a  combined  output  of  213  dB  rel  1  pPa@lm. 

A  variety  of  57.0  Hz  continuous  wave  (CW)  and  phase  encoded  signals  were 
transmitted  on  a  fixed  schedule  consisting  of  one  hour  of  continuous  transmission 
followed  by  two  hours  of  silence  repeated  for  the  entire  test  period.  The  transmission 
schedule  is  included  in  Appendix  A.  A  carrier  frequency  of  57.0  Hz  was  chosen  to  avoid 
excessive  attenuation  at  higher  frequencies  and  increased  shipping  noise  at  lower 
frequencies.  In  addition,  the  resonance  frequency  of  the  transducers  was  close  to  57.0  Hz. 
Figure  2.1  shows  the  relation  between  predicted  ambient  noise  spectra  and  the  transmitter 
frequency  response. 

The  three  types  of  transmitted  signals,  described  in  detail  in  Chapter  HI,  are 

as  follows: 

•  Continuous  Wave  at  57.0  Hz. 
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•  3/10  Pentaline.  This  is  a  phase  encoded  signal  with  five  major  spectral  lines  of 
various  amplitude  which  facilitates  easy  estimates  of  signal-to-noise  ratio  based 
upon  the  observed  frequency  lines.  However  only  minimal  line  resolution  is 
possible  since  the  power  is  split  evenly  between  the  signal  and  the  carrier. 

•  M-Sequences.  Maximal  length  shift  register  sequences,  sometimes  called  pseudo¬ 
random  noise,  are  phase  encoded  signals  of  various  lengths  which  correspond  to  a 
specific  code.  Like  pentaline,  only  half  the  power  is  in  die  signal.  However,  pulse 
compression  techniques  are  used  to  provide  excellent  time  resolution. 


Hz 


Figure  2.1:  Average  Ambient  Noise  Level  and  Transmitter  Frequency  Response 
[Ref.  8]. 
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3.  Reception 

Since  the  vertical  array  in  Bermuda  suffered  a  disabling  casualty,  the  array 
deployed  jointly  by  the  Naval  Postgraduate  School  (NPS),  Monterey  Bay  Aquarium 
Research  Institute  (MBARI),  Massachusetts  Institute  of  Technology  (MIT),  and  Woods 
Hole  Oceanographic  Institute  (WHOI)  served  as  the  only  participating  vertical  array  and 
as  such  will  be  the  sole  domain  of  this  thesis.  Figure  2.2  illustrates  the  optimal 
positioning  of  the  vertical  array,  based  on  ray  tracing  (Ort  and  Chiu,  1990)  utilizing  a 
global  ocean  mesoscale  variability  model  (Semtner  and  Chervin,  19S8).  This  computer 
simulation  demonstrated  that  acoustic  rays  emanating  from  the  Heard  Island  transmission 
source  at  angles  greater  than  127°  true  would  pass  to  the  east  of  New  Zealand  and  that 
rays  less  than  135°  would  escape  interaction  with  Antarctica.  In  addition,  these  rays 


13 


had  five  or  less  interactions  with  the  ocean  bottom.  This  band  of  rays,  shown  in  Figure 
2.2,  intercepts  the  west  coast  of  the  United  States  in  the  Monterey  coastal  region  between 
San  Francisco  and  Point  Conception,  making  it  the  most  likely  area  on  the  U.S.  west 
coast  to  receive  the  signal  because  of  the  convergence  of  rays  in  this  vicinity. 

The  32  hydrophone  array  was  suspended  vertically  from  R/V  Point  Sur  at 
approximately  36°  30’  N  and  122°  45’  E,  as  shown  in  Figure  2.3.  The  hydrophones 
were  spaced  45  meters  apart  with  the  top  one  nominally  positioned  at  a  depth  of  345 
meters  and  the  bottom  one  at  1740  meters.  3000  meters  of  kevlar  core  cable  was  used 
to  house  50  pair  of  #28  solid  copper  conductors.  The  inboard  100’  of  cable  was  shielded 
with  copper  to  minimize  induced  currents  from  the  R/V  Point  Sur.  The  lower  2  km  had 
SAIC  quiet  cable  coating  over  a  28"  twist  reversal  pattern  of  0.8”  per  side  which  reduced 
strumming  of  the  suspended  cable.  The  600  pound  (submerged  weight)  array  and  its 
suspended  ballast  weight  of  300  lbs.  were  balanced  by  the  buoyancy  from  a  main 
subsurface  float  and  approximately  half  of  the  30  "football"  floats,  spaced  2  meters  apart, 
which  isolated  the  array  from  the  surface  gravity  wave  field.  The  main  subsurface  float 
was  a  37"  sphere  with  475  lbs.  of  buoyancy.  The  football  floats  measured  7.5"  x  15"  and 
had  a  buoyancy  of  6.25  lbs.  each.  The  lower  weight  was  a  55  gallon  mass  attached  with 
an  acoustic  release  to  facilitate  release  at  recovery.  The  midpoint  of  the  array  was 
positioned  with  about  neutral  buoyancy  near  the  center  of  the  deep  sound  channel  at  about 
1100  meters.  An  instrument  to  record  PSK  telemetry,  pressure,  array  tilt,  temperature, 
and  conductivity  was  attached  approximately  4  meters  above  the  top  hydrophone  An 
internally  recording  pressure  and  array  tilt  sensor  was  attached  just  below  the  lowest 
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hydrophone.  The  submerged  portion  of  the  array  was  distanced  from  the  vessel  by  about 
1500  meters  using  4  large  surface  floats  spaced  by  20  smaller  floats. 


B.  ACOUSTIC  SIGNAL  PROPAGATION 

1.  Propagation  Theory 

Sound  propagation  in  a  fluid  can  be  described  by  the  linearized,  lossless  wave 
equation  (Kinsler  and  others,  1982).  In  a  motionless  medium,  the  equation  is 


v?p  =  -L  Up 

c2  8  t2 


(2.4) 


where  c(x,y,z)  is  the  sound  phase  speed  field  and  p(x,y,z)  is  the  fluid  pressure  disturbance 
field.  The  linearized,  lossless  wave  equation  is  a  good  approximation  in  a  motionless 
medium.  Variations  in  the  travel  time  due  to  current  are  generally  an  order  of  magnitude 
less  than  the  changes  due  to  deviations  in  the  sound  speed  and  so  can  be  ignored  in  a  first 
order  treatment  (Spindel,  1986).  As  an  example,  consider  a  mesoscale  eddy  in  which  a 
typical  10  cm/s  surface  current  diminishes  to  about  3  cm/s  below  the  thermocline.  A 
corresponding  temperature  anomaly  near  the  sound  channel  axis  is  approximately  1  degree 
Celsius.  This  temperature  perturbation  results  in  a  sound  speed  perturbation  of 
approximately  5  m/s  which  is  100  times  the  change  due  to  current.  Also  note  that  a 
change  in  the  path  length  (which  could  occur  due  to  mesoscale  variability)  can  have  an 
even  greater  effect  on  the  travel  time,  depending  on  the  size  of  the  change  (Flatte  and 
others,  1979). 
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In  discussing  sound  wave  propagation,  it  is  useful  to  define  a  propagation 

vector 


1c  =  kx  x  +  kyy  +  kzz 

that  has  a  magnitude  (or  wave  number)  given  by 


(2.5) 


that  gives  the  location  of  the  point  (x,y,z)  with  respect  to  the  origin  v0,0,0)  of  a 
rectangular  coordinate  system.  We  can  let 


*  =  C°*  -  =  tl  <*•*> 

w  K0 

where  c0  is  a  constant  reference  speed  and  *¥  has  units  of  length.  Note  that  if  c  =  Cj,  then 
'F  =  (k,x  +  k,y  +  kxz)  /  K  which  is  the  distance  from  the  origin.  Surfaces  of  constant 
'F  will  now  correspond  to  surfaces  of  constant  phase.  That  is,  the  locus  of  all  points 
(x,y,z)  satisfying 

vjf  =  constant  +  c0 1  (2.9) 

defines  a  surface  of  constant  phase  at  each  instant  t,  as  it  propagates  through  the  acoustic 

medium.  VT  is  always  perpendicular  to  surfaces  of  constant  phase  and  is  therefore  the 

local  direction  of  propagation  of  the  phase  surface.  Substitution  into  the  wave  equation 
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with  A(x,y,z)  representing  the  pressure  amplitude  field  yields 


V2v4 

— )vt|rVi|i  +  (— T 

~j 

— +  V2t|») 

A 

(co)  UJ 

col  A  ) 

=  0 


(2.10) 


When  the  real  part  (the  value  within  the  first  square  brackets)  is  set  equal  to 
zero,  the  result  is  the  Eikonal  equation.  When  the  imaginary  part  (the  value  within  the 
second  square  brackets)  is  set  equal  to  zero,  the  Transport  equation  results.  If  A  and  V'P 
satisfy  the  inequalities 


V2A 

A 

UJ 

|V2  i|r|  «  —  ,  and 


VM 


Vt|f 


<  «  (2.11) 

c 


then  all  the  terms  in  Eq.  2.10  except  the  second  and  third  can  be  neglected  and  the 
Ei1'onal  equation  can  be  approximated  by 


Vty  -  V  t(r  = 


<c  '2 
V  c  / 


=  n2 


(2.12) 


where  n(x,y,z)  is  the  refractive  index  as  a  function  of  position.  This  is  known  as  the 
plane  wave  approximation  of  the  wave  equation.  The  inequalities  of  Eq.  2.11  can  be 
described  as  (Kinsler  and  others,  1982) 

1 .  The  amplitude  of  the  wave  must  not  change  appreciably  in  distances  comparable 
to  the  wavelength. 

2.  The  speed  of  sound  must  not  change  appreciably  in  distances  comparable  to  the 
wavelength. 

3.  The  channel  thickness  and  source-receiver  distance  must  be  large  in  comparison 
to  a  wavelength. 
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The  57.0  Hz  signal  used  in  this  experiment  had  a  wavelength  of  about  26.3  meters  which 
satisfies  all  the  above  conditions  and  can  thus  be  described  using  the  plane  wave 
approximation.  A  solution  to  the  plane  wave  equation  can  be  expressed  as 

p  =Ae^ot~*'r)  (2.13) 

The  essence  of  ray  theory  is  that  V'F  defines  the  direction  of  wave 

propagation  as  a  function  of  position.  Therefore  a  point  by  point  solution  of  V*P  renders 
the  path  traversed  by  the  acoustic  energy.  This  technique  has  been  implemented  with  the 
Hamiltonian  raytracing  program  HARPO  (Jones  and  others,  1986)  and  interfaced  with 
output  from  a  global  eddy  resolving  ocean  general  circulation  model  (Semtner  and 
Chervin,  1990),  to  predict  the  optimum  receiver  locations  for  the  Heard  Island  Feasibility 
Experiment  (Chiu  and  Ort,  1990).  Typical  three  dimensional  results  are  illustrated 
horizontally  in  Figure  1.4  and  vertically  in  Figure  2.4.  The  rays  emanating  from  the 
Heard  Island  source  at  true  horizontal  angles  between  127°  and  135°  were  seen  to  avoid 
significant  interaction  with  the  bottom,  land  masses,  and  seamounts  that  extend  into  the 
sound  channel.  These  propagating  rays  can  be  divided  into  two  general  classifications 
depending  on  their  behavior  in  the  vertical  plane. 

Type  1  rays,  shown  in  Figure  2.4(a),  are  trapped  close  to  the  surface  by  a 
shallow  surface  duct  until  the  duct  weakens  at  ranges  between  7500  and  8000  km.  The 
rays  then  enter  the  deep  sound  channel.  Because  the  distance  between  the  surface  duct 
and  the  deep  sound  channel  is  almost  1200  m,  the  rays  enter  the  deep  sound  channel  at 
a  relatively  steep  angle  and  then  oscillate  about  the  axis  with  a  large  spatial  amplitude. 
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Type  2  rays  are  also  trapped  near  the  surface  in  the  half  channel  regime,  but 
less  so  than  the  Type  1  rays.  At  a  range  of  about  5000  km,  the  Type  2  rays  encounter 
a  stronger  vertical  gradient  which  refracts  the  acoustic  energy  downward  toward  the  deep 
sound  channel.  Also  at  about  5000  km,  the  deep  sound  channel  axis  gradually  increases 
depth  to  about  1000  m  over  a  distance  of  about  600  km.  This  gradual  entry  angle  into 
the  deep  sound  channel  causes  the  Type  2  rays  to  have  a  much  smaller  vertical 
fluctuations  about  the  axis. 

Since  the  Type  1  rays  oscillate  further  from  the  minimum  sound  speed  at  the 
channel  axis,  they  traverse  water  with  a  higher  sound  speed  during  most  of  their  path. 
They  can  thus  be  expected  to  arrive  sooner  than  the  axially  trapped  Type  2  rays.  In 
addition,  note  that  the  upper  portion  of  the  Type  1  ray  trajectory  lies  within  the  oceanic 
thermocline  where  sound  speed  variability  is  the  greatest.  Therefore,  although  the  Type 
1  rays  are  more  representative  of  the  upper  ocean  (which  is  more  closely  coupled  to  the 
atmospheric  temperature),  they  will  also  exhibit  more  variability  in  arrival  time. 

The  travel  time  for  a  specific  ray  X\  (the  1th  ray  path)  can  be  found  by 
integrating  the  sound  slowness  (inverse  speed)  over  the  path 

t,  =  f  -  (2.14) 

Jx,  c 

The  perturbations  in  sound  speed  can  then  be  viewed  as  deviations  from  some  arbitrary 
reference  speed  c„, 

c  =  Cq  +  6  c  ,  1 2.15) 

so  that  the  total  travel  time  becomes  a  constant  travel  time  with  a  perturbation 
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(2.16) 


H  (2.18) 

co 

2.  Propagation  Losses 

Little  prior  experience  exists  with  propagation  losses  over  the  18  Megameter 
ranges  traversed  in  this  experiment.  However  expanding  the  current  theories  to  global 
ranges  can  lend  some  insight  to  the  order  of  magnitude  losses  that  can  be  expected.  The 
most  significant  energy  loss  should  be  spreading  losses  due  the  expanding  area  of  the 
phase  surface  as  the  wave  propagates  from  the  source.  A  secondary  loss  due  to 
absorption  and  scattering  can  also  be  expected. 


a.  Spreading  Losses 

A  worst  case  scenario  involves  totally  lossy  boundaries  which  can  be 
approximated  by  spherical  spreading  in  an  iso-velocity  fluid  (Birdsall,  1990).  Invoking 
energy  conservation  in  a  lossless  medium  on  an  acoustic  phase  surface  propagating 
outward  from  a  point  source  gives  the  spherical  wave  solution  at  range  r  from  the  sv»urce 
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(2.19) 


p  =  -^2—  eJ<.o<  -  *0 
4  nr2 

It  follows  that  the  energy  contained  on  a  sphere  of  unit  diameter  about  the  point  source, 
p0,  is  conserved  on  the  expanding  phase  surface  as  shown  by 

p0  -  4nr12/]  =  4  nr22/2  (2.20) 

where  I,  and  I2  are  the  phase  surface  intensities  at  ranges  1  and  2  respectively.  The 

transmission  loss  is  then  given  by 

TL  =  10  log—  =  101ogr22  =  201ogr2  (2.21) 

4 

which  yields  a  loss  of  145.1  dB  over  18  Mm. 

A  best  case  scenario  assumes  totally  reflective  boundaries  which  can  be 
approximated  by  a  flat  earth  similar  to  a  right  cylinder  with  an  18  Mm  radius  and  5  km 
height  (Birdsall,  1990).  The  cylindrical  wave  solution  is 

n  =  en»>  -  in  (2.22) 

4nr 

The  transmission  loss  is  then  given  by  10  times  the  log  of  the  intensity  on  the  cylindrical 
phase  surface  relative  to  a  unit  sphere 

(2nr2H)  ( r,H ) 

TL  =  lOlog - —  =  lOlod—  (2.23) 

V".’j 

where  r,  is  1  m  and  H  is  the  depth  of  water  in  meters.  This  formula  yields  a 
transmission  loss  of  106.5  dB  over  18  Mm.  This  flat  earth  model  assumed  2n  x  18  Mm 
=  113.1  Mm  as  the  cylindrical  circumference  although  the  circumference  is  physically 
constrained  to  the  40  Mm  circumference  of  the  earth.  Birdsall  suggests  using  a 
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convergence  factor  of  10  log(l  13.1/40)  =  4.5  dB  to  correct  for  this  difference  which 
reduces  the  predicted  transmission  loss  to  102  dB.  If  the  depth  is  constrained  to  a 
nominal  1000  m  height  of  the  sound  channel,  the  predicted  transmission  loss  is  reduced 
to  99.5  dB.  Both  models  are  pessimistic  in  that  some  of  the  observed  intensities  will  be 
greater  than  the  average  predicted  intensity.  However,  both  models  are  also  optimistic 
in  that  they  assume  no  interaction  with  the  surface,  the  bottom,  or  with  oceanic  fronts. 

b.  Absorption  Losses 

Absorption  at  57  Hz  is  usually  considered  negligible,  but  over  18  Mm  the 
cumulative  loss  becomes  significant.  The  57  Hz  carrier  frequency  falls  below  the  peak 
relaxation  frequencies  of  both  MgS04  and  B(OH)3  but  is  dominated  by  the  B(OH)3 
relaxation  process  (Clay  and  Medwin,  1977).  For  a  depth  of  about  1000  m  and  a 
temperature  of  4°  C,  we  can  find  an  absorption  coefficient  a  from  (Thorp,  1967) 

a  -  -  *  0.000275/2  +  0.0007  (2-24) 

1  *f2  4.1 

where  f  is  the  carrier  frequency  and  a  is  in  dB  per  kyd.  The  terms  in  the  above  equation 
are  due  to  boric  acid  ionization,  magnesium  sulphate  relaxation,  viscosity  shear,  and  low 
frequency  scattering  respectively.  The  last  term  is  dominant  at  57  Hz  and  is  taken  from 
Figure  2.5  which  compares  the  empirical  values  of  low-frequency  attenuation  to  the 
extrapolated  relaxation  absorption.  An  absorption  coefficient  of  0.0007  dB  per  kyd  or 
0.758  dB  per  Mm  gives  an  absorption  loss  of  13.6  dB  over  18  Mm. 
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Frequency,  kHz 


Figure  2.5:  Low  Frequency  Absorption  Coefficient  [Ref.  15}. 

3.  Signal-to  Noise  Ratio 

The  signal  level  at  the  receiver  hydrophones  can  be  predicted  by 

Received  Signal  Level  *  Source  Level  -  Transmission  Loss  (2.25) 
where  all  levels  are  interpreted  as  band  levels.  This  gives  levels  between  54.3  anc  97.4 

dB  rel  1  p Pa  at  the  receiver  hydrophones  depending  on  the  spreading  model  used. 
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The  predicted  57  Hz  noise  at  the  receiver,  shown  in  Figure  2.1,  is  primarily 
due  to  shipping  and  falls  between  75  to  93  dB  rel  1  pPa.  The  signal-to-noise  level  can 
be  predicted  by 

SNR  =  Receiver  Signal  Level  -  Receiver  Noise  Level  (2.26) 

which  gives  an  SNR  from  -38.07  to  22.4  dB  depending  on  the  assumed  scenario.  Using 

these  values,  between  0  and  56.93  dB  of  array  and/or  processing  gain  is  needed  to  meet 
our  arrival  time  resolution  criteria  of  18.86  dB.  The  method  of  maximizing  the 
processing  gain  will  be  described  in  the  following  section. 

C.  MAXIMAL  LENGTH  SEQUENCES 

1.  Introduction 

Impulse  excitation  (a  signal  of  infinitesimal  length  and  infinite  magnitude)  is 
a  convenient  and  powerful  method  for  measuring  the  response  of  a  system  and  for 
determining  the  travel  time  through  a  media.  The  oceanic  acoustic  medium  can  be  treated 
as  a  large,  time-varying  distortionless  filter  whose  response  to  impulsive  excitation  is 
simply  the  sum  of  delayed  impulses  traveling  along  the  individual  paths  (Dees,  1989) 

Paths 

hit )  =  J2  ax  hit  -  x)  (2-27) 

<=! 

where  a,  is  the  amplitude  of  the  i0*  path,  and  x,  is  the  travel  time  along  the  i"*  path. 

Unfortunately,  ideal  impulse  transmissions  are  physically  difficult  to  achieve. 
The  amplitude  and  length  of  a  pulse  are  limited  by  the  peak  power  and  bandwidth, 
respectively,  of  the  transmitter.  Impulsive  acoustic  signals  can  be  generated  by  explosive 
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or  implosive  sources  but  these  have  an  uneven  frequency  distribution  of  energy,  and  are 
generally  not  reproducible.  Another  solution  is  to  use  pseudorandom  noise  to  provide  a 
means  of  pulse  compression.  In  this  method,  a  long  coded  signal  is  transmitted  whose 
period,  frequency  distribution,  and  energy  are  deterministic  and  can  be  tailored  to  meet 
experimental  requirements.  The  received  signal  can  then  be  cross-correlated,  through  a 
matched  filter,  with  the  original  code  to  produce  a  single  pulse  with  a  much  higher 
amplitude  and  shorter  length  than  the  transmitted  signal.  In  addition,  the  coded  periodic 
signal  can  be  repeated  for  additional  processing  gain. 

The  pulse  compression  technique  used  for  the  Heard  Island  Feasibility 
Experiment  was  to  utilize  a  phase-modulated  carrier  signal  to  transmit  a  specific  sequence 
of  ones  and  zeros  known  as  a  maximal  length  sequence  or  m-sequence  (Ziemer  and 
Peterson,  1985).  The  most  important  characteristic  of  the  m-sequence  is  its 
autocorrelation,  which  is  constant  except  at  a  shift  of  zero  (corresponding  to  perfect 
correlation  with  the  transmitted  code),  making  the  sequence  analogous  to  white  noise. 
The  width  of  the  correlated  pulse  is  equal  to  the  width  of  one  individual  digit  of  the  code, 
Tc,  which  was  0.087  seconds  for  this  experiment.  Details  of  the  specific  signal  used  in 
this  experiment  are  given  in  Chapter  m. 

The  amplitude  at  the  correlated  zero-shift  digit  has  an  increase  over  amplitude 
of  the  coded  signal  equal  to  the  number  of  digits  in  the  code,  N.  The  ideal  pulse 
compression  gain  for  a  single  sequence  of  code  is  given  by  10  log(N).  Processing  gain 
is  thus  dependent  on  the  code  length  which  is  in  turn  limited  by  the  length  of  tin  e  the 
ocean  can  be  assumed  stationary.  One  of  the  goals  of  this  experiment  was  to  determine 
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the  maximum  effective  sequence  length.  This  experiment  used  sequence  lengths  of  255, 
511,  1023,  and  2047  digits,  corresponding  to  22.37,  44.8,  89.7,  and  179.6  seconds 
respectively.  The  ideal  processing  gain  for  these  sequences  was  then  24.06,  27.08,  30.1, 
and  33.11  dB  respectively.  In  addition,  if  the  m-sequence  code  is  repeated  L  times,  a 
gain  of  10  log(L)  can  be  added  to  the  single  sequence  gain  under  ideal  conditions.  The 
amount  of  correlated  repetition  is  dependent  on  the  stability  of  the  arrival  structure  which 
is  in  turn  influenced  by  the  meso,  gyre,  and  basin  scale  variability  of  the  ocean 
environment. 

2.  M-Sequence  Generation 

M-Sequences  are  generated  using  a  simple  binary  shift  register  whose 
structure  is  defined  by  a  specified  polynomial  (Ziemer  and  Peterson,  1985).  A  general 
shift  register  generator  is  shown  in  Figure  2.6.  A  deterministic  and  periodic  sequence 
of  bits  can  be  output  from  any  one  of  the  registers,  with  period  NTC.  Although  an  equally 
valid  m-sequence  is  generated  in  both  the  forward  and  reverse  order,  this  experiment 
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consistently  used  the  sequential  output  from  the  0th  register,  R*,  such  that  the  consecutive 
bits  generated  at  the  output  comprised  the  m-sequence.  In  addition,  although  the  initial 
state  of  the  shift  register  is  arbitrary,  this  experiment  consistently  used  an  initial  sequence 
of  all  zeroes  followed  by  a  one. 

The  shift  register  structure  is  specified  by  a  primitive  polynomial  such  as 

gD  -  bfi’  *  bfi  *  f>0  «•»> 

where  n  is  the  number  of  delays  in  (or  degree  of)  the  shift  register  generator,  D  is  the  unit 

delay  and  b„  are  the  feedback  weighting  coefficients.  The  weighting  coefficients  take  on 

values  of  one  or  zero  indicating  connection  or  no  connection,  respectively,  to  the  n,h 

register.  Table  8-5  from  Ziemer  and  Peterson  lists  octal  representations  (referred  to  as 

the  octal  law)  of  the  binary  coefficients  corresponding  to  particular  generating 

polynomials.  M-sequences  do  not  repeat  a  register  state  until  after  2"-l  delays.  The 

governing  polynomial  is  known  as  a  primitive  polynomial  since,  for  a  specified  number 

of  shift  register  stages,  a  sequence  of  maximal-length  period  is  produced  before  the 

sequence  repeats.  The  signal  specifications  for  the  m-sequence  code  used  in  this 

experiment  are  detailed  in  Chapter  HI.  In  practice,  the  m-sequence  digits  are  transformed 

by  mapping  [1,0]  to  [-1,1]  which  enables  use  of  faster  correlation  procedures.  Because 

many  people  are  familiar  with  binary  mathematics,  however,  all  arithmetic  operations  for 

polynomials  in  this  thesis  are  performed  using  finite-field  arithmetic  (modulo  two). 

3.  Characteristics 

The  output  of  an  n-stage  m-sequence  shift  register  generator  with  digit  length 
Tc  and  length  N  =  2°  - 1  can  be  represented  as 


f 
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(2.29) 


C  =  ai  ^t-i Te 

where  a,  =  (-l)1*,  and  U  is  a  unit  pulse  beginning  at  0  and  ending  at  Tc.  Its  discrete  two- 
sided  autocorrelation  function  Rk  for  k  cyclic  shifts  of  the  generator  is  given  by 


1.0  for  k-mN\ 


K  =  l 


(2.30) 


—  for  k*mN 
N 


where  m  is  any  integer. 


The  non-discrete  autocorrelation  function  of  the  m-sequence  code  for  time 


delay  of  t  is  given  by 


T,  T, 

R  =  1  -  -3  Rk  +  —  *A+1 
\  c)  \  c) 
where  xt  -  x  -  (N  -  1)TC  (Ziemer  and  Peterson,  1985) 

If  0  <  t  <  Tc  ,  then  k  =  0  and  xt  =  x  so  that 


(231) 


R  -  1  -  — 

c  T. 


1  T 


N  T. 


-  4  ♦ 

tI  n) 


(232) 


If  Tc  £  x  <,  (N-l)Te ,  then  k  *  mN  and  (k+1)  *  mN  for  any  integer  m  so 


T  \  1  1 

R  =  1  -  — 

c  T  N  T  N 


(233) 


If  (N-1)TC  <  x  <  (NTJ  ,  then  k  =  N  -  1  and  k  +  1  =  1  so  that 
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(2.34) 


Now  Eqs.  32,  33,  and  34  define  one  complete  cycle  of  the  maximal  length  sequence 
autocorrelation  which  is  illustrated  in  Figure  2.7. 


R 


Figure  2.7:  Autocorrelation  function  for  a  m-sequence  of  length  N. 

The  Wiener-Khinchine  theorem  shows  that  the  Fourier  transform  of  an 
autocorrelation  function  describes  the  frequency  distribution  of  power  and  is  referred  to 
as  the  power  spectral  density  or  power  spectrum  of  the  time  series  (Bendat  and  Piersol, 
1986).  The  power  spectrum  of  the  baseband  m-sequence  is  then 

S,  -  r  A,  i(f  -  if„)  <2-35) 

where  A,  =  1/N2,  A,  =  [(N+l)/N2]  sinc2(i/N),  and  fN  =  1/NTC  This  spectrum  appears  as 
a  series  of  discrete  impulse  lines  separated  by  the  code  repetition  frequency  1/NT(  with 
the  line  amplitude  envelope  defined  by 
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(236) 


A  =  «^2(/tc) 

l  N1) 

for  all  except  the  DC  line  which  has  an  amplitude  of  1/N2. 

If  the  phase  shift  angle  0  is  chosen  such  that 

0  =  tan_1(v^V)  (237) 

then  the  spectrum  of  the  carrier  signal  will  fall  exactly  on  the  envelope  of  impulses  and 

result  in  the  maximum  signal -to-noise  ratio  after  demodulation  and  sequence  removal 
(Spindel,  1985).  The  optimum  phase  shifts  for  pulse  compression  are  then  approximately 
86.42°,  87.47°,  88.21°,  and  88.73°  for  the  255,  511,  1023,  and  2047  digit  sequences 
respectively.  The  transmitted  signal  for  this  experiment  used  a  90°  phase  change  between 
digits  for  all  modulations  which  provided  roughly  half  the  power  in  the  carrier.  This 
phase  shift  was  chosen  because  it  produces  about  6  dB  steps  in  the  power  spectrum  of 
the  pentaline  signal  which  facilitates  graphical  estimation  of  the  signal-to-noise  ratio  as 
discussed  in  the  following  Chapter. 

When  the  m-sequence  c  is  used  to  biphase  modulate  a  carrier  wave  having 
power  Pc  and  frequency  f0  the  resultant  signal  is  given  by 

$  =  y/2  Pc  c  cos(2nf0t)  (238) 

The  power  spectrum  of  the  modulated  carrier  signal  S,  is  then  the  convolution  of  the 

carrier  power  spectrum  Sc  and  the  m-sequence  power  spectrum 

S,  *  ♦  f  «</  -  /»>  *  k  *  \  »</ +  /„) 

and  the  effective  result  is  a  translation  of  the  discrete  spectrum  Sf  upward  and  downward 
by  a  frequency  f0.  The  spectrum  of  the  transmitted  signal  is  then  a  band  of  impulses 
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centered  at  f0  and  separated  by  1/NTC.  Figure  2.8  illustrates  the  amplitude  envelope  of 
the  individual  impulses  from  the  transmitted  signal. 
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Figure  2.8:  Envelope  of  Power  Spectrum  impulses  in  the  Transmitted  M-Sequence 
Signal. 


4.  Autocorrelation  Process 

Maximal  length  sequences  have  several  properties  which  make  it  attractive  as 
a  pulse  compression  code.  Its  autocorrelation  function  has  a  large  magnitude  and  short 
duration  which  enhances  arrival  time  estimates.  It  is  also  deterministic  in  that  once  the 
polynomial  is  specified,  all  output  parameters  are  known.  Additionally,  it  produces  a 
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repeatable,  periodic  signal  which  can  be  easily  implemented.  The  single  major  drawback 
is  that  standard  autocorrelation  procedures  are  computationally  intensive. 

The  sample  autocorrelation  function  of  a  stationary  (ergodic)  random 
process  with  a  zero  mean  can  be  estimated  as  (Spindel,  1985) 


*zx  =  dt  OiX  it  (2*40) 

1  T  0 

for  time  t  and  time  delay  t.  For  N  data  values  {x.},  n=l,2,...N,  sampled  at  equally  spaced 
time  intervals  5t  from  a  stationary  record  x^  with  a  zero  mean,  R„  can  be  numerically 
computed  by 


1 

n  -  r 


N-r 


E 


r  =  0,l,2,...,/n 


(2.41) 


where  F  is  the  lag  number  and  m  is  the  maximum  lag  number  (m  <  N).  Since  the  arrival 
time  is  not  known  in  advance,  standard  correlation  methods  require  the  input  signal  to  be 
correlated  with  all  N  shifted  versions  of  the  original  sequence.  Thus  N2  multiplications 
are  necessary  for  computing  a  standard  FFT  based  cross-correlation.  The  following 
paragraphs  describe  a  less  computationally  intensive  autocorrelation  method  which  takes 
advantage  of  the  m-sequence  structure. 

Using  a  3  register  generator  output  as  an  example  with  an  octal  law  of  13  and 
a  polynomial  of  degree  3,  the  forward  signal  code  S  is 


5  =  [1001110]  <2-42> 

To  correlate  S  with  all  seven  shifted  versions  of  itself,  we  can  construct  a  matrix  M  of 

all  possible  versions 
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M  = 


10  0  1110 
0  10  0  111 
10  10  0  11 
110  10  0  1 
1110  10  0 
0  1110  10 
0  0  1110  1 


(2.43) 


When  the  digits  of  S  and  M  are  mapped  [1,0]  to  [-1,1],  multiplying  the  signal  by  the 
correlation  matrix  will  result  in  the  autocorrelation, 

Rm  =  MS  (2.44) 

This  product  is  the  entire  goal  of  the  initial  signal  processing  but  requires  a  computer 
intensive  N2  add-multiply  operations.  A  faster,  more  efficient  algorithm,  called  the  Fast 
Hadamard  Transform  can  be  used  to  reduce  the  number  of  operations  required  to  perform 
the  autocorrelation. 

5.  The  Hadamard  Matrix 

To  facilitate  discussion  of  the  Fast  Hadamard  Transform,  it  is  necessary  to 
introduce  the  Hadamard  matrix  (Ziemer  and  Peterson,  1985).  The  Sylvester-type 
Hadamard  matrix  has  a  recursive  form  for  higher  orders  given  by 


1  1 

Ht  Ht 

Hx  =  [1]  , 

sc: 

II 

1  "I 

7 

Hi+\ 

For  our  third  order  example,  the  matrix  H  is 
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It  is  useful  to  map  the  Hadamard  matrix  from  [1,-1]  to  [0,1]  as  follows: 

0  0  0  0  0  0  0  0 
0  10  10  10  1 
0  0  1  1  0  0  1  1 
0  110  0  110 
H  ~  0  0  0  0  1  1  1  1 
0  10  110  10 
0  0  11110  0 
0  110  10  0  1 


(2.46) 


(2.47) 


H  can  be  factored  into  a  binary  counting  (from  1  to  7  in  this  example)  matrix  A  and  its 
transpose  AT,  which  is  also  a  binary  counting  matrix  as  given  by 
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H =  AAT  = 


0  0  0 
0  0  1 
0  1  0 
0  1  1 
1  0  0 
1  0  1 
1  1  0 
1  1  1 


0  0  0  0  1  1  1  1 
0  0  1  1  0  0  1  1 
0  10  10  10  1 


(2.48) 


The  correlation  matrix  M  can  be  similarly  factored  into  two  matrices  B  and 
C  which  contain  all  binary  combinations,  but  not  as  simply.  The  first  matrix  B  is 
formed  from  the  sequential  output  of  the  shift  register,  but  bit  reversed  (from  right  to  left) 
and  in  the  reverse  order  (from  bottom  to  top).  The  original  order  is  then  preserved  by 
shifting  the  rows  of  the  matrix  to  bring  the  3x3  identity  matrix  to  the  top.  For  our  third 
order  m-sequence  example,  the  factor  matrix  B  is  given  by 


B  = 


1  0  0 
0  1  0 
0  0  1 
1  1  0 
0  1  1 
1  1  1 
1  0  1 


(2.49) 


The  second  matrix  C  is  formed  from  any  3  shifted  versions  of  the  m-sequence,  ie 


C  = 


10  0  1110 
0  10  0  111 
10  10  0  11 


(2.50) 
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It  is  easily  verified  that 


BC  =  M  (2*51) 

By  adding  a  leading  column  of  zeros  to  B,  forming  B\  and  adding  a  leading  row  of  zeros 
to  C,  forming  C’,  a  new  matrix  M’  can  be  formed  by 

M'  =  B'  C'  ^2*52) 

Note  that  M*  contains  all  possible  combinations  of  ones  and  zeros,  as  does  A,  albeit  in 

a  different  order.  If  mapping  matrices  can  be  found  such  that  QA  =  B’  and  ATP  =  C\ 

then  the  same  matrices  will  map  the  Hadamard  matrix  to  the  m-sequence  matrix 

M'  =  B'C'  =  QAATP  =  QHP  ^2-53) 

6.  The  Fast  Hadamard  Transform 

Recall  that  correlation  of  the  received  signal  S  with  the  m-sequence  code  M 
was  performed  by  the  matrix  product  of  S  and  M.  By  adding  a  leading  zero  to  S, 
forming  S’,  equation  (2.44)  becomes 

R  sm  =  M'S>  *  QHPS'  (2-$4) 

Equation  (2.S4)  shows  that  matrix  M  is  not  needed  and  the  Hadamard  matrix  can  be  used 

in  its  place.  An  efficient  method  of  forming  the  matrix  product  of  QHPS’  is  the  Fast 

Hadamard  Transform  (FHT).  The  FHT  takes  advantage  of  the  fact  that  any  vector  PS’ 

which  is  multiplied  by  the  unmapped  Hadamard  matrix  (consisting  of  +l’s  and  -l’s), 

results  in  a  sum  of  all  the  terms  of  PS’  with  various  +  and  -  weighting.  This  summation 

can  be  effectively  utilized  to  cancel  random  noise  amplitudes  where  a  random  sum  of 

positive  and  negative  values  will  equal  0.  In  addition,  non-random  amplitudes  contained 

in  the  m-sequence  code  will  sum  coherently  to  a  value  equal  to  the  sequence  length  times 
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the  absolute  amplitude  of  a  single  value.  For  our  third  degree  example  an  arbitrary  vector 
given  by 


a 

b 


c 


e 


f 

8 

h 


(2.55) 


after  multiplication  by  the  Hadamard  matrix  becomes 


'1 

1 

1 

1 

1 

1 

1 

>1 

a 

a  +6 +c  +d+e  +f+g+h 

1  -1 

1  -1 

1  -1 

1  -1 

b 

a-b+c-d+e-f+g-h 

1 

1 

-1  -1 

1 

1 

-1  -1 

c 

a+b-c-d+e+f-g-h 

HPS'  = 

1  -1 

-1 

1 

1  -1 

-1 

1 

d 

m 

a-b-c+d+e-f-g+h 

1 

1 

1 

1 

-1  -1 

-1  -1 

e 

a+b+c+d-e-f-g-h 

1  -1 

1  -1 

-1 

1 

-1 

1 

f 

a-b+c-d-e+f-g+h 

1 

1 

-1  -1 

-1  -1 

1 

1 

8 

a+b-c-d-e-f+g+h 

i  -• 

-1 

1 

-1 

1 

1  -1 

h 

a  -b -c  +d-e  +f+g -h 

(2.56) 


which,  when  permuted  by  Q,  is  the  autocorrelation  function  R,,.  This  technique  is  called 
the  Fast  Hadamard  Transform.  When  written  in  the  form  of  a  flow  graph  as  in  Figure 
2.9,  the  product  HPS’  has  the  same  form  as  the  well  known  fast  Fourier  transform,  where 
the  complex  multiplications  are  set  to  one  (Cohn  and  others,  1976  and  Borish  and  others, 
1983).  Note  that  there  is  no  bit  reversal  or  multiplication  by  a  phase  factor  as  is  required 
by  the  Fast  Fourier  Transform.  Since  the  multiplication  by  P  and  Q  has  been  replaced 
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Basic  FHT  Element 


a+b+c+d+e+f+g+h 
a-b+c-d  +  e-f+g-h 
a+b-c-d+e+f-g-h 
a-b-c+d+e-f-g+h 
a+b+c+d-e-f-g-h 
a-b+  c-d-e+  f-g+  h 
a+b-c-d-e-f+g+h 
a-b-c+d-e+f+g-h 


Figure  2.9:  The  Basic  Fast  Hadamard  Transform  Element  for  Cascading 
Additions  and  the  Full  Diagram  for  an  8  Point  FHT. 


by  permuting  the  matrix  positions,  there  is  no  multiplication  required.  The  increase  in 
speed  of  execution  is  difficult  to  compare  with  the  FFT  since  this  method  only  requires 
additions.  During  the  correlation  procedure,  a=0  so  no  new  information  is  added.  The 
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zeroth  position  result  can  be  seen  to  be  the  sum  of  all  the  elements  of  the  code  and  is 
thus  equal  to  the  DC  pedestal.  This  pedestal  can  be  removed  by  subtracting  the  zeroth 
position  from  all  other  positions.  All  that  remains  is  to  determine  the  permutation 
matrices  P  and  Q. 

7.  Vector  Order  Permutation 

The  permutation  matrices  P  and  Q  must  be  found  such  that  QA  =  B’  and  ATP 
=  C\  Matrices  A,  AT,  B’,  and  C’  contain  all  possible  arrangements  of  ones  and  zeros 
in  3  digit  rows  or  columns.  A  natural  index  for  each  row  or  column  is  then  its  equivalent 
octal  value.  The  indices  are  important  in  that  they  allow  ’multiplication’  by  the 
permutation  matrices  by  rearranging  the  order  of  the  signal  vector,  rather  than  direct 
multiplication.  Note  that  no  multiplications  are  required,  only  the  reordering.  For  a  given 
code,  the  permutations  can  be  evaluated  once,  prior  to  signal  processing,  and  the  result 
stored  as  an  index  array  to  be  applied  to  each  vector. 

8.  Correlation  Summary 

The  following  is  a  summary  of  the  autocorrelation  procedure: 

1.  Augment  the  signal  vector  S  with  a  zero  in  the  zeroth  position  to  form  S’. 

2.  Permute  the  augmented  input  vector  S’  according  to  P. 

3.  Perform  the  Fast  Hadamard  Transform  additions  as  indicated  by  the  flow  graph 
in  Figure  2.9. 

4.  Permute  the  vector  resulting  from  the  FHT  according  to  Q. 

5.  Remove  the  zeroth  entry  and  subtract  it  from  all  other  elements  to  remove  the  DC 
bias. 
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III.  SIGNAL  PROCESSING 

A.  SIGNAL  SPECIFICATION 

1.  Continuous  Wave 

The  CW  signal  was  a  pure  sinusoidal  waveform  transmitted  at  a  nominal  array 
power  level  of  213  dB  rel  luPa@lm,  This  signal  was  useful  in  estimating  the  frequency 
stability  of  the  arrivals,  and  in  comparing  signal  levels  to  noise  in  the  frequency  band  of 
interest.  However,  the  extreme  power  levels  at  which  the  transducers  were  driven  coupled 
with  severe  weather  conditions  in  the  vicinity  of  Heard  Island  caused  variability  in  the 
output  level  performance  of  the  transmitter. 

2.  3/10  Pentaline 

The  3/10  pentaline  signal  has  five  major  spectral  lines  hence  the  name 
pentaline.  It  consists  of  3  digits  of  10  carrier  cycles  each  at  57.0  Hz,  a  total  of  30  cycles 
per  period,  repeated  continuously  for  one  hour.  This  produces  a  nominal  bandwidth  of 
57.0  Hz  /  10  cycles  per  digit  =  5.7  Hz.  The  signal  is  specified  by 

s  =  A  cos  (2n fct  +  bnQ)  ,  (3*1) 

where  b„  =  II:  -1,  +1,  +1  :  II,  corresponding  to  law  7octtJ  or  digits  of  1,0,0  in  binary  and  n 
is  incremented  once  every  10  carrier  cycles.  The  first  digit  is  phase  shifted  by  -45 
degrees,  the  second  and  third  are  phase  shifted  by  +45  degrees.  There  are  57.0  H;:  /  30 
cycles  per  period  =  1.9  periods  per  second,  which  produces  a  spectral  line  separation  of 
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Frequency  (hz) 


Figure  3.1:  Power  Spectrum  of  the  Pentaiine  Signal. 

1.9  Hz.  Figure  3.1  illustrates  the  power  distribution  in  the  frequency  domain.  The  power 


levels  rel  57.0  Hz  CW  signal  are  tabulated  below. 


LOWER  FREQ. 

|  (Hz) 

UPPER  FREQ. 
(Hz) 

POWER 
(Percent  of  CW) 

POWER 
(dB  rel  CW) 

57.0 

57.0 

55.5 

-2.5 

55.1 

58.9 

14.7 

-8.3 

53.2 

60.8 

3.5 

-14.5 

51.3 

62.7 

0.0 

-OO 

49.4 

64.6 

0.8 

-20.8 

47.5 

66.5 

0.5 

-229 
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Note  that  the  power  decreases  in  approximately  6  dB  steps  in  line  amplitude  which 
permits  graphical  estimation  of  received  signal-to-noise  ratios  from  the  observed  power 
spectrum. 


3.  Maximal  Length  Sequences 

The  maximal  length  sequences  each  consist  of  a  series  of  digits  of  3  cycles 
each  at  57.0  Hz,  producing  a  digit  duration  of  87.7  milliseconds  and  a  bandwidth  of  11.4 
Hz.  The  signal  is  specified  as 

s  =  Acos(2n/cf  +  bnQ)  ,  (3*2) 

where  n  is  incremented  once  every  5  carrier  cycles.  Phase  modulation  is  +/-  45  degrees 
as  it  is  for  the  3/10  pentaline.  The  initial  digits  are  0,0, ...0,0,1  i.e.,  the  maximum  string 
of  0’s  followed  by  a  1.  The  sequence  lengths  are  255,  511,  1023,  or  2047  digits  with 
their  characteristics  tabulated  below. 


Digits 

LAW 

(octal) 

PERIOD 

(seconds) 

Periods/Hr 

255 

537 

22.36 

160.9 

511 

1473 

44.82 

80.3 

1023 

2033 

89.82 

40.0 

2047 

5747 

179.56 

MNHKSftl 

B.  DATA  COLLECTION  AND  PRE-PROCESSING 

1.  Amplification 

The  receiving  array  of  32  hydrophones  continuously  collected  32  channels  of 
data  measuring  the  acoustic  field  at  depths  between  approximately  345  and  1740  meters. 
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The  hydrophones  had  a  nominal  sensitivity  of  -170  dB  (rel  lV/pPa)  corresponding  to  316 
pPa/pV.  Therefore,  the  data  can  easily  be  converted  directly  to  pressure  expressed  in 
micropascals.  Received  hydrophone  data  underwent  a  linear  fixed  gain  of  20  dB  at  a 
patch  panel  and  then  was  linearly  preamplified  for  a  gain  of  6  dB  before  processing 
through  a  variable  gain  amplifier.  The  variable  gain  amplifier  automatically  switched  gain 
(1,  8,  84,  512,  or  4096)  depending  on  the  input  level  to  optimize  the  10  volt  dynamic 
range  of  the  analog  to  digital  converter  and  thus  the  precision  of  the  recorded  signal  data. 
All  data  analyzed  in  this  thesis  has  been  normalized  to  the  hydrophone  output  voltage 
levels  using  the  program  convert. c  (listed  in  Appendix  B)  which  is  an  implementation  of 
the  program  noimize.c  written  by  Von  der  Heydt  of  Woods  Hole  Oceanographic 
Institution  (per  E-mail  communication  dtd  17  April  1991). 

2.  Analog  Filtering 

As  noted  above,  the  maximum  bandwidth  for  the  transmissions  is  1 1.4  Hz  for 
all  of  the  m-sequence  signals.  It  is  desirable  to  bandpass  filter  the  received  data  to 
eliminate  any  out-of-band  noise.  Therefore,  in  addition  to  amplification,  the  signal 
underwent  bandpass  filtering  prior  to  recording.  The  highpass  filter  had  a  -6  dB  point  at 
10  Hz  with  a  rolloff  of  -30  dB/octave.  The  lowpass  filter  had  a  -3  dB  point  at  80  Hz 
with  a  rolloff  of  -48  dB/octave. 

3.  Analog-to-Digital  Conversion 

Figure  3.2  illustrates  the  data  collection  and  pre-processing  apparatus  al>oard 
R/V  Point  Sur.  The  acoustic  signal  was  digitally  sampled  at  228  Hz  for  transmissions  of 
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Figure  3.2  Data  Preprocessing 


26  January  1991  to  January  1991  and  then  at  456  Hz  from  January  to  February  1991. 
These  sampling  frequencies  are  4  and  8  times  the  Nyquist  frequency  respectively  and  are 
an  integer  multiple  of  the  optimum  carrier  frequency  fc.  Sampling  in  this  manner 
simplifies  processing  and  interleaving  of  data,  since  it  produces  an  integer  number  of  data 
points  per  digit.  The  data  was  digitally  recorded  continuously  during  the  1  hour  reception 
intervals  onto  both  1  Gigabyte  optical-disk  and  VHS  tape  media.  After  normalization  to 
hydrophone  output  levels,  the  data  was  stored  in  an  array  of  binary  single  precision 
floating  point  values.  Real  time  processing  and  signal  monitoring  was  conducted  using 
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both  a  Concurrent  6400  Real  Time  Unix  processor  and  more  conventional  spectrum 
analyzers. 

4.  Clipping 

During  postprocessing  of  the  data,  occasional  broadband  spikes  of  high 
amplitude  noise  of  undetermined  origin  were  noticed  which  lasted  approximately  one 
second.  These  noise  transients  tended  to  saturate  the  local  time  series  to  a  degree  that 
attempts  to  recover  the  signal  near  the  spike  were  futile.  Since  these  spikes  in  the  time 
domain  data  had  no  relation  to  the  transmitted  signal  but  significantly  raised  the 
background  noise  level,  they  can  be  clipped  without  significantly  degrading  the  received 
signal.  During  preprocessing,  if  the  absolute  value  of  a  data  point  exceeded  twice  the 
standard  deviation  from  the  mean  of  its  local  10  second  buffer  it  was  clipped  or  limited 
to  that  value.  By  clipping  to  a  standard  deviation  from  a  localized  mean,  short  term  noise 
anomalies  are  eliminated  but  the  general  characteristics  of  signal-to-noise  in  the  slowly 
varying  ambient  noise  field  are  retained.  Clipping  was  incorporated  into  the  resampling 
module  but  occurred  prior  to  resampling. 

5.  Doppler  Correction 

Both  the  transmitting  and  receiving  ships  and  arrays  were  acted  upon  by  the 
forces  of  tides,  currents,  winds,  etc.,  which  combined  to  produce  some  amount  of  relative 
drift  between  them.  In  addition,  the  source  array  was  towed  during  transmission  periods 
as  required  to  maintain  steerage  of  the  source  vessel  with  minimum  thrust.  Although  both 
vessels  were  operated  to  minimize  array  motion  (especially  accelerations),  some  relative 
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movement  was  unavoidable  and  because  of  this  the  received  signal  center  frequency  was 
often  Doppler  shifted  away  from  the  transmitted  center  frequency.  Since  the  motions 
were  slow  and  generally  uniform,  it  is  sufficient  to  make  corrections  for  first  order 
Doppler  shifts  and  neglect  acceleration  effects  (Birdsall  and  others,  1987,  1988). 

When  the  signal  s  incurs  a  Doppler  frequency  shift  due  to  relative  velocity, 
v,  the  received  signal  experiences  an  expansion  or  compression  in  time.  The  Doppler 
shifted  signal  sd  can  be  expressed  as 

5<* =  s[(i  *  9*  -  t<]  (33) 

where  c  is  the  average  phase  speed  of  sound  in  water  given  in  meters  per  second,  and  t, 
is  the  delay  of  the  i**  arrival  from  the  transmitter  to  the  receiver.  Since  the  objective  of 
the  Heard  Island  Feasibility  Experiment  was  to  examine  the  resolution  of  the  arrival 
structure  and  not  to  determine  the  travel  time  itself,  T,  is  arbitrary  and  only  important  in 
its  relation  to  other  arrival  delays.  Therefore  it  will  be  assigned  a  reference  value  of  zero 
for  this  discussion.  The  magnitude  of  the  Doppler  shifted  power  spectrum  is  given  by 


(3.4) 


When  the  Doppler  shifted  signal  is  digitally  sampled  at  a  frequency  of  f,  the  value  of  the 
n,b  sample  is  given  by 


(3.5) 


If  the  sampling  frequency  f,  is  adjusted  by  the  same  Doppler  shift  a  new  sampling 
frequency  results  and  is  given  by 
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(3.6) 


which  would  produce  a  new  n"1  sample 

’  S(±J}' 

When  Eq.  3.6  is  substituted  into  Eq.  3.7,  the  resampled  signal  is 


(3.7) 


(3.8) 


which  is  identical  to  a  signal  with  no  Doppler  shift  and  thus  signal  processing  can 
proceed  as  in  the  zero  Doppler  case. 

Prior  to  signal  processing,  each  data  set  was  resampled  at  either  4  or  8  times 
the  Doppler  shifted  carrier  frequency  for  the  data  originally  sampled  at  228  Hz  or  456  Hz 
respectively.  The  re-sampled  values  resulted  from  linear  interpolation  between  actual  data 
points  as  if  the  desired  sampling  frequency  had  been  used.  The  average  received  center 
frequency  was  estimated  by  averaging  repeated  32,768  point  (140.35  second)  power 
spectrums  of  the  data  over  the  one  hour  transmission  period.  The  resultant  time-averaged 
peak  frequency  nearest  57.0  Hz  was  used  as  the  received  center  frequency.  This  value 
agrees  with  subsequent  estimations  of  Doppler  shift  using  the  drift  of  the  source  and 
receiver  ships  from  GPS  data  and  the  relative  geodesic  angle.  Since  changes  in  Doppler 
were  small,  no  correction  was  made  for  variations  in  received  frequency  during  the 
reception  hour. 
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C.  SIGNAL  PROCESSING  PROCEDURE 

1.  General  Considerations 

Data  processing  was  performed  using  a  proprietary  signal  processing 
applications  software  package  named  Lab  Workbench  which  was  developed  by 
Concurrent  Computer  Corporation.  This  software  provides  a  shell  interface  between  a  32 
channel  analog-to-digital  converter,  the  Concurrent  6400  Real-Time  Unix  based  computer 
system,  numerous  signal  processing  software  modules,  and  the  optical  disk  data  storage 
medium.  Most  importandy,  Lab  Workbench  allows  optional  user-defined  data  processing 
modules  which  are,  in  effect,  Fortran  or  C  language  subroutines  which  can  operate  real¬ 
time  on  the  data  as  it  is  being  collected. 

Because  hardware  difficulties  prevented  the  use  of  the  real-time  capabilities 
during  the  Heard  Island  Feasibility  Experiment,  all  data  was  processed  after  the 
conclusion  of  the  experiment,  using  a  playback  module  to  read  the  data  from  the  storage 
media.  The  stored  data  was  retrieved  from  the  optical  disk  and  processed  in  buffer 
lengths  of  1  second  each.  The  1  second  buffers  of  data  were  passed  sequentially  through 
a  bandpass  filter  module,  a  demodulator  module,  two  parallel  low-pass  filter  modules,  a 
sequence  removal  module,  and  finally  a  module  which  writes  the  processed  data  onto  the 
storage  media.  Each  of  the  signal  processing  modules  as  well  as  pre-processing  steps  are 
described  in  sequential  order  below.  Figure  3.3  illustrates  the  time  domain  processing 
accomplished  with  the  Lab  Workbench  modules.  Listings  of  C  language  source  code  for 
the  processor  modules  are  included  in  Appendix  B.  Modular  design  of  the  processing 
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Figure  3.3  Overview  of  Signal  Processing  Steps 
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software  allowed  increased  flexibility  of  implementation  and  operation  in  any  one  of  the 
following  modes: 

•  Setup:  Specification  of  all  run-time  variables,  memory  allocation,  and  preliminary 
calculations.  Performing  routine  calculations  and  ‘housekeeping’  functions  in  this 
mode  speeds  operations  during  the  run  mode. 

•  Run:  Actual  real-time  data  processing  which  occurs  in  buffer  lengths  of  one  second 
each  until  the  entire  one  hour  transmission  has  been  processed. 

•  Stop:  Halts  the  Run  phase  of  data  processing  but  maintains  current  data  in  the 
transfer  buffers  to  allow  subsequent  restart. 

•  Menu:  Allows  selection  of  variable  values  such  as  carrier  frequency,  m-sequence 
code,  phase  shift,  etc.  during  operation  or  execution  of  the  Lab  Workbench 
software. 

•  Drop:  Same  function  as  Stop  but  drops  the  last  data  transfer  buffer. 


2.  Bandpass  Filter  Module 

The  power  spectrum  of  the  phase  encoded  m-sequence  is  the  convolution  of 
the  carrier  power  spectrum  and  the  m-sequence  code  power  spectrum  as  developed  in 
Section  2.C.3.  Therefore  the  retrieved  data,  which  was  hardware  bandpass  filtered  30 
to  80  Hz  prior  to  recording  was  again  bandpass  filtered.  This  second  bandpass  filter 
passed  the  band  from  45  to  69  Hz  using  a  software  implementation  of  a  classical 
Hamming  windowed,  20th  order,  linear-phase,  finite  impulse  response  (FIR)  filter. 
Received  data  was  passed  in  a  24  Hz  band  corresponding  to  the  power  spectrum  of  the 
biphase  modulated  57  Hz  carrier.  The  filter  structure  is  the  general  tapped  delay-line 
filter  described  by  the  difference  equation 
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(3.9) 


y„  =  b0Xn  +  bxxn~l  +  •  ■  •  +b20Xn-20 

for  data  vector  x  and  weighting  coefficients  b,  obtained  using  the  firl  function  in  Matlab. 
Matlab  is  a  proprietary  software  product  developed  by  The  MathWoiks,  Inc.,  South 
Natick,  MA.  Although  the  C-language  Lab  Woricbench  implementation  of  this  filter  is 
listed  in  Appendix  B  as  bandpass.c,  any  other  filter  type  can  logically  be  substituted  if 
it  maintains  linearity  in  phase. 

3.  Demodulation 

To  retrieve  the  digital  code  from  the  biphase  modulated  carrier  signal,  the 
received  data  must  be  demodulated.  The  object  of  complex  or  quadrature  demodulation 
is  to  isolate  the  real  bandpass  signal  on  the  positive  frequency  axis  and  translate  it  to  a 
complex  lowpass  signal.  If  fd  is  the  downshift  frequency  in  Hz  and  t  is  time  in  seconds, 
the  demodulation  can  be  accomplished  by  multiplying  the  input  data  by  a  complex  rotator 
given  by  (Birdsall  and  Metzger,  1990) 

d,  =  ste~i2*fdt  (3*10) 

where  s,  represents  the  input  signal  and  dt  is  the  demodulated  output  signal.  When  the 

sampling  rate  f,  is  precisely  the  fourth  harmonic  of  the  carrier  frequency  fe,  as  it  is  for  the 

re-sampled  data,  fdt  changes  by  exactly  one-quarter  cycle  with  each  successive  sample  and 

the  values  of  the  downshifting  complex  rotator  e'i2*<fc/f,)n  being  l,-i,-l,i;l,  -i,  -1,  i;  etc.  The 

complex  multiplication  then  results  in  successive  polarity  changing  and  alternating 

assignment  to  real  and  imaginary  parts.  The  result  is  a  real  and  imaginary  train  of  pulses 

corresponding  to  the  m-sequence  code  at  a  baseband  of  0  to  11.4  Hz.  Since  the 

demodulator  is  not  synchronized  with  the  transmitter,  it  must  demodulate  the  fignal 
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without  the  benefit  of  phase  information.  Therefore,  to  recover  all  of  the  magnitude  of 
the  signal  in  the  baseband,  the  signal  must  be  multiplied  by  both  the  cosine  and  sine  at 
the  carrier  frequency.  The  digital  implementation  of  the  real  part  of  the  demodulation  is 
given  by 


r  =  s.cos 

n  n 


(2  nfen 


(3.11) 

l  f.  I 

where  r.  is  the  n*  demodulated  real  data  point.  A  listing  of  the  C-language  code  is  given 
in  Appendix  B  as  demod.c.  The  digital  formulation  of  the  imaginary  part,  i„,  is  similarly 


iB  =  snsm\ 


(3.12) 

Ideally,  the  downshifting  frequency  fd  would  be  the  received  signal  carrier  frequency 
which  would  translate  the  signal  spectrum  perfectly  to  it’s  baseband.  However,  since  the 
signal  was  re-sampled  at  4  or  8  times  the  measured  fe,  we  can  use  any  combination  of  fc 
and  f,  in  equations  3.11  and  3.12  such  that  f/fc=4  or  8  and  still  maintain  correct 
synchronization  for  the  complex  demodulation  rotator.  Since  the  transmitter  design  carrier 
frequency  was  57.0  Hz,  the  ideal  sampling  frequency  of  either  228  Hz  or  456  Hz  were 
used  for  fc  and  f,  respectively. 


4.  Low  Pass  Filter 

Both  the  in-phase  and  quadrature  components  from  the  demodulator  were  low 
pass  filtered  passing  the  band  from  0  to  12  Hz.  The  filters  were  a  software 
implementation  of  a  classical  Hamming  windowed,  20th  order,  linear-phase,  finite  impulse 
response  (FIR)  filter.  The  passed  band  from  0  to  12  Hz.  corresponded  to  the  power 
spectrum  of  the  baseband  of  the  biphase  modulated  57  Hz  carrier.  The  filter  structure  is 
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again  the  general  tapped  delay-line  niter  described  by  the  difference  equation  (Gq.  3.9) 
with  vector  x  and  coefficients  b„  obtained  using  the  firl  function  in  Matlab.  Although 
the  software  implementation  of  this  filter  is  listed  in  Appendix  B  as  lowpass.c,  any  other 
filter  type  can  logically  be  substituted  if  it  maintains  linearity  in  phase. 

5.  Sequence  Removal 

Correlation  of  the  received  data  with  the  transmitted  maximal  length  sequence 
code  is  performed  by  a  software  adaptation  of  the  Fast  Hadamaid  Transform.  The 
permutation  matrices  discussed  in  Chapter  II  are  calculated  for  the  specific  m-sequence 
code  and  stored  during  the  setup  mode  of  the  sequence  removal  module. 

During  run-time  execution  of  the  module,  dual  inputs  of  real  and  imaginary 
data  are  input  and  multiplexed  into  a  large  storage  array.  Since  the  Hadamaid  processing 
must  be  performed  on  an  entire  sequence,  the  storage  array  will  hold  22.4,  44.8,  89.8,  or 
179.6  seconds  of  data  prior  to  processing,  depending  on  the  specific  m-sequence  code. 
Remember  that  the  transmitted  signal  had  5  cycles  per  digit  and  the  signal  was  sampled 
at  4fc  or  K  so  there  were  20  data  points  collected  per  digit  for  the  228  Hz  sampled  data 
and  40  points  per  digit  for  the  456  Hz  sampled  data.  When  the  signal  was  demodulated 
the  number  of  data  points  per  digit  grew  by  a  factor  of  two,  since  the  demodulated  signal 
has  real  and  imaginary  parts.  The  incoming  real  r  and  imaginary  i  data  buffers  are 

interleaved  into  the  storage  array  by  alternating  r0,  i„,  r„  i„ . r39n,  i39n  for  an  n-digit 

sequence  sampled  at  228  Hz. 

Since  the  Fast  Hadamard  Transform  only  requires  one  data  point  per  digit,  by 
the  time  a  single  m-sequence  has  been  received  the  storage  array  actually  contains  to  or 
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80  interleaved  m-sequences  depending  on  the  sampling  rate.  Each  sequence  is 
sequentially  permuted  according  to  the  scrambling  array,  processed  by  the  FHT,  and  then 
permuted  again  according  to  the  unscrambling  array.  When  the  Fast  Hadamard  Transform 
is  performed  successively  on  each  interleaved  sequence  in  the  storage  array,  the  result 
from  an  ideal  signal  is  40  or  80  adjacent  peaks  at  die  0-shift  or  correlated  digit.  This 
oversampling  causes  a  flat-topped,  one  digit  wide  correlation  peak  vice  the  ‘pointy’ 
autocorrelation  peak  developed  in  Chapter  II. 

As  shown  in  Figure  2.7,  the  autocorrelation  function  for  an  m-sequence  has 
a  negative  DC  level  for  all  digit  positions  except  the  first.  Therefore,  the  zeroth  position 
of  each  sequence  contains  the  average  EXT  level  and  can  be  calculated  and  removed  from 
all  other  positions.  The  complex  level  adjustment  Cw  is  related  to  the  average  DC  level 

L««  by 

-  /,  <3J3> 

and  f,  is  given  by  the  relation 

f  =  j  (W+l)tan(6)  +  j(N-tan2(Q))  (3 14) 

1  N2  +  tan2(0) 

where  N  is  the  m-sequence  digit  length  and  0  is  the  phase  modulation  angle.  All  the 
parameters  of  equation  3.15  are  known  and  so  f,  is  calculated  only  once  during  the  setup 
phase  of  the  sequence  removal  module.  Bias  removal  then  requires  one  complex 
multiplication  per  sequence  during  the  run  mode  to  calculate  C^  and  the  correctin'1  can 
be  subtracted  from  each  value  in  the  sequence. 
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Since  the  processed  data  array  is  interleaved  with  real  and  imaginary  parts  of 
the  transformed  signal,  it  is  a  trivial  matter  to  compute  the  magnitude  mag.  in  the  time 
domain  by 

MSn  =  fiT+V  (3*15) 

The  ideal  processed  signal  now  is  a  time  series  which  has  zero  values  at  all  points  except 
the  correlated  digit.  At  the  correlated  digit,  there  will  be  a  sequence  of  pulses  whose 
combined  width  is  one  digit  wide  (0.087  seconds)  and  whose  amplitudes  are  equal  to  the 
number  of  digits  in  the  sequence.  For  real  data,  the  model  output  is  similar  with  the 
addition  of  randomly  distributed  white  noise  present  in  all  digits  as  shown  in  Figure  3.2. 

Although  a  Doppler  shifted  input  signal  can  be  correctly  demodulated  without 
resampling,  the  bandpass  digit  sequence  would  still  be  either  expanded  or  compressed  in 
time.  The  Fast  Hadamard  Transform  scheme  used  in  this  thesis  assumes  exactly  20  digits 
per  sample  (40  for  the  456  Hz  sampled  data)  and  only  one  sample  per  m -sequence  digit 
is  transformed  at  a  time.  Thus  any  expansion  or  contraction  of  the  time  series  greater 
than  one  sample  interval  would  cause  a  loss  of  a  digit.  A  single  output  of  the  FHT 
process  would  have  a  reduction  in  amplitude  of  1/N  for  each  shift  of  one  sample  length 
where  N  is  the  number  of  digits  in  the  sequence.  Because  the  interleaved  data  produces 
20  (or  40)  adjacent  output  peaks,  the  reduction  occurs  first  on  the  end  peaks  and  then 
successively  toward  the  opposite  end  as  the  Doppler  shift  increases.  The  result  is  a 
narrowing  of  the  pulse  at  the  top  with  increasing  Doppler  and  an  eventual  decrease  in  the 
pulse  amplitude  for  Doppler  time  shifts  exceeding  one  digit  width  Te  of  0.0877  seconds 
over  the  length  of  the  sequence.  The  Doppler  frequency  shifts  required  to  cause  the  loss 
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of  one  complete  digit  range  from  0.028  Hz  (1.43  kts)  for  the  2047  sequence  to  0.2224  Hz 
(1 1.5  kts)  for  the  255  sequence.  Therefore,  it  is  imperative  that  the  data  is  resampled  to 
remove  Doppler  time  alteration  and  preserve  a  constant  number  of  samples  per  digit. 

D.  POST-PROCESSING  ENHANCEMENTS 

1.  Time  Averaging 

The  primary  objective  of  time  domain  processing  is  to  obtain  a  correlated 
pulse  whose  amplitude  is  sufficiently  greater  than  the  average  noise  level  to  allow 
confident  detection.  Since  the  m-sequence  is  transmitted  repeatedly  during  the  one-hour 
transmission  period,  the  repetitive  sequences  can  be  averaged  incoherently  by  simply 
summing  the  corresponding  data  points  of  successive  sequences.  As  the  number  of 
averages  increases,  a  stable  correlated  signal  pulse  will  remain  at  the  same  relative 
position  in  the  time  series  but  the  Guassian  noise  will  tend  toward  an  average  value  of 
zero.  The  sequence  removal  module  facilitates  averaging  by  retaining  leftover  buffer  data 
when  a  full  sequence  has  been  received  so  that  it  can  be  used  at  the  beginning  of  the 
following  sequence.  In  this  way  a  continuous  stream  of  data  is  processed  as  it  was 
received. 

2.  Visual  Integration 

Because  the  received  signal  is  often  very  weak  even  after  significant 
processing  to  maximize  gain,  it  is  useful  to  display  the  resultant  data  in  a  manner 
conducive  to  easy  interpretation  and  recognition  of  the  signal.  It  is  particularly  helpful 
to  display  the  graphical  output  of  several  time  series  such  that  corresponding  times  are 
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aligned.  In  this  way,  any  signal  which  may  be  marginally  detectable  in  a  single  time 
series  may  stand  out  as  a  trend  among  several  sequential  series. 

The  Fortran  program  hafton.f  which  is  listed  in  Appendix  B  was  written  to 
take  advantage  of  this  visual  integration  between  sequential  time  series.  It  produces  a 
two-dimensional  plot  with  m-sequence  time  series  on  the  horizontal  axis  and  successive 
m-sequences  ‘stacked’  above  in  sequential  order.  The  amplitudes  are  represented  by 
halftones  of  various  gray  shades  with  higher  amplitudes  indicated  by  darker  shading. 
These  plots  can  be  further  enhanced  by  plotting  a  function  of  the  data  such  as  log,  sine, 
cosine,  etc.  to  accentuate  the  extreme  or  mean  regions  of  data  amplitudes. 

An  additional  method  of  achieving  visual  integration  is  the  use  of  a  ‘waterfall’ 
display  wiiere  successive  time  or  frequency  domain  plots  are  ‘stacked’  and  recurring 
peaks  become  apparent.  The  matlab  program  mesh.m  was  used  to  generate  the  waterfall 
type  displays  included  in  this  thesis. 
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IV.  RESULTS 


A.  DATA  COLLECTION 

The  first  objective  of  this  thesis  was  to  successfully  deploy  a  vertical  array  in  the 
optimum  reception  zone  and  record  the  sampled  acoustic  environment  at  a  rate  higher 
than  the  Nyquist  frequency.  This  goal  was  met  cooperatively  between  the  Naval 
Postgraduate  School,  the  Monterey  Bay  Aquarium  Institute,  Woods  Hole  Oceanographic 
Institution,  and  Massachusetts  Institute  of  Technology.  The  output  of  the  transducers  was 
sampled  at  either  228  or  456  Hz,  bandpass  filtered,  and  optimally  amplified  to  maximally 
utilize  the  dynamic  range  of  the  analog-to-digital  converter  as  discussed  in  Chapter  HI. 

At  01261525  (1525  Greenwich  Mean  Time  on  26  January  1991)  R/V  Point  Sur  was 
on  station  near  the  northern  edge  of  the  predicted  reception  zone,  recording  all 
transmissions  until  the  array  cable  became  entangled  in  the  ships  screw  at  01300205 
GMT.  During  this  period  the  R/V  Point  Sur  recorded  12  transmissions  and  steamed 
between  reception  periods  from  35°  41’N,  123°  01’W,  to  35°  30’N,  123°05’W  in  order 
to  sample  a  cross  section  of  the  predicted  reception  zone.  Following  a  swift  transit  to 
Monterey  in  order  to  free  the  fouled  screw  and  an  admirable  cable  splicing  repair  effort 
led  by  K.  Von  der  Heydt  of  WHOI  and  K.  Lashkari  or  MBARI,  R/V  Point  Sur  returned 
to  a  point  further  into  the  predicted  reception  zone  at  34°  15’N,  122°  00’W  by  01312108 
GMT.  Unfortunately,  by  this  time  the  transmitting  vessel,  RV  Cory  Chouest,  had  become 
encumbered  by  heavy  seas  and  the  reliability  of  it’s  overdriven  transducers  was 
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deteriorating  rapidly.  Two  abbreviated  receptions  were  recorded  between  01312108  GMT 
and  02010039  GMT  at  the  new  location  before  untenable  conditions  near  Heard  Island 
forced  conclusion  of  the  Experiment.  The  32  channels  of  digitized  data  were  sampled  and 
recorded  to  1  Gigabyte  Tahiti  ZCAV  optical  disks  at  58.386  kilobytes  per  second  for  the 
228  Hz  sample  rate  and  1 16.736  kilobytes  per  second  for  the  456  Hz  sample  rate. 

Daily  sound  speed  measurements  were  recorded  using  an  XBT  and  representative 
sound  velocity  profiles  are  shown  in  Figure  4.1.  In  addition  array  depth  and  tilt  data  was 


Sound  Speed  (m/s)  Sl*cd 


Figure  4.1:  Sound  Velocity  Profiles. 
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recorded  from  upper  and  lower  instruments  on  the  array.  This  data  is  used  to  determine 
hydrophone  location  relative  to  the  sqund  channel  axis  located  at  the  sound  velocity 
minima.  The  top  hydrophone  was  nominally  at  345  meters  depth.  Follow-on  data 
analysis  will  use  the  array  depth  and  tilt  data  to  implement  modal  beamforming  with  all 
32  hydrophones.  Beamforming  can  theoretically  provide  10  log(32)  or  15  dB  of  array 
gain  and  a  realistic  expectation  of  between  6  and  10  dB.  In  addition,  if  the  hydrophone 
outputs  are  properly  weighted,  the  array  will  act  as  a  mode  filter  to  discriminate  between 
various  normal  modes  of  propagation  and  more  accurately  depict  the  propagation  method. 

B.  SIGNAL  PROCESSING 

1.  Objective 

The  second  goal  of  this  thesis  was  to  develop  and  implement  signal  processing 
software  capable  of  optimizing  the  characteristics  of  the  transmitted  signal.  The  software 
was  designed  to  operate  in  real  time  on  an  incoming  data  stream  using  a  Concurrent  6400 
data  acquisition  computer  with  Lab  Workbench  applications  software.  To  this  end,  the 
programs  in  Appendix  B  which  were  discussed  in  Chapter  HI  have  been  developed.  All 
processing  was  implemented  with  user  function  modules  in  the  Lab  Workbench 
environment.  All  user  function  modules  were  written  in  the  C  programming  language  and 
can  be  utilized  real-time  as  the  data  is  being  collected.  The  32  channel  input  data  stream 
can  be  demultiplexed  with  a  channel  selector  switch  to  permit  analysis  of  individual 
channels  or,  if  desired,  arithmetically  combined  to  sum  individual  channels.  Piping  data 
between  modules  is  facilitated  by  connecting  the  input  and  output  boxes  of  appropriate 
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modules  on-screen  with  a  click  of  the  left  mouse  button.  In  addition,  user  interface  with 


the  running  program  was  enhanced  with  on-screen  menus  to  permit  selection  of  key 
parameters  for  each  function.  The  on-screen  menus  are  activated  by  clicking  the  left 
mouse  button  while  the  cursor  is  positioned  on  the  module  of  interest. 

2.  Test  Data 

The  signal  processing  programs  were  tested  with  an  artificial  data  set  which 
was  recorded  onto  audio  cassette  tape  and  then  played  back  through  a  single  channel  of 
the  data  collection  system.  Test  data  consisted  of  66  consecutive  sequences  of  255  digit 
m-sequence  code  (octal  law  537).  The  inherent  noise  of  the  audio  tape  and  variability  of 
the  tape  drive  speed  provided  both  a  source  of  near-white  noise  and  frequency  wavering 
similar  to  the  Doppler  shift  of  the  drifting  source  and  receiver.  The  results  of  the 
performance,  and  sensitivity  testing  using  the  test  data  are  demonstrated  below. 

3.  Resampling  and  Clipping 

Figure  4.2  compares  the  signal  band  power  spectrum  of  the  01270005  m- 
sequence  raw  data  to  the  same  data  after  resampling  and  clipping.  The  C  function 
resamp.c  clips  and  resamples  32  channels  of  multiplexed  data  at  a  rate  exactly  four  times 
the  Doppler  shifted  carrier  frequency  as  discussed  in  Chapter  HI.  The  function  clips  the 
data  at  a  magnitude  which  is  a  user  defined  factor  times  the  standard  deviation  (a)  of  a 
local  10  second  time  series.  It  allows  the  user  to  define  both  the  Doppler  shifted 
resampling  frequency  and  the  clipping  factor  from  an  on-screen  menu.  All 
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Frequency  (Hz) 


7  xlQ"  (fa-1  RESAMPLED  AND  CLIPPED  POWER  SPECTRUM  OF  01270003  CH-6 


Frequency  (Hz) 


Figure  4.2:  Results  of  Clipping  and  Resampling  01270005  Raw  Data. 


data  analyzed  for  this  thesis  was  clipped  at  the  default  value  of  2 o  and  resampled  at 
exactly  4  times  the  median  carrier  frequency  as  determined  by  spectral  analysis.  A 
32,768  point  FFT  was  used  in  all  spectral  analysis  in  this  thesis.  Resampling  simplifies 
implementation  of  the  demodulation  procedure  by  minimizing  the  effects  of  relative 
motion  induced  Doppler  frequency  shifts  of  the  received  carrier  frequency.  The  clipping 
of  this  particular  time  series  resulted  in  approximately  a  3  dB  decrease  in  the  median 
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4.  Bandpass  Filtering 


Figure  4.3  displays  the  filter  response  of  the  20th  order  FIR  bandpass  filter  for 
both  228  and  456  Hz  sample  rates  superimposed  on  the  normalized  magnitude  envelope 
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Figure  4.3:  Magnitude  and  Phase  Response  of  20th  Order  Bandpass  Filter, 
of  the  m-sequence  signal.  The  filter  passes  frequency  components  between  45  and  79  Hz. 
It  is  readily  apparent  that  removal  of  unwanted  out  of  band  noise  will  significantly 
improve  the  signal-to-noise  ratio  of  the  processed  data  since  unwanted  noise  in  the  time 
domain  is  comprised  of  a  wide  spectrum  of  frequencies.  Bandpass  filtering  was 
performed  with  the  function  bandpass.c.  The  program  allows  user  selection  of  the  sample 
frequency  (which  determines  the  filter  coefficients)  from  an  on-screen  menu. 

5.  Demodulation 

The  complex  demodulation  of  the  received  signal  into  its  quadrature  and 
cophase  components  at  baseband  frequency  was  also  implemented  with  a  Lab  Workbench 
user  function  module.  It  is  listed  in  Appendix  B  as  fdemodl.c.  Both  the  real  and 
imaginary  parts  of  the  signal  are  produced  as  output  streams  for  further  manipulation, 


65 


thereby  preserving  all  of  the  signal  information  for  maximum  processing  gain.  The  user 
can  select  carrier  frequency  and  sample  frequency  from  an  on-screen  menu.  Figure  4.4 
shows  a  1  second  portion  of  the  test  data  time  series  after  processing  with  the 
demodulator.  This  one  second  series  contains  11.4  digits  each  of  length  0.0877  seconds. 


Time  (seconds) 


Figure  4.4:  Real  Part  of  Demodulated  255  M -Sequence  Test  Data. 

6.  Lowpass  Filtering 

The  output  of  the  demodulator  in  Figure  4.4  exhibits  some  high  frequency  (57 
Hz)  ripple  components  in  the  time  series.  Since  the  m-sequence  has  been  translated  to 
its  baseband  frequency  (0  to  11.4  Hz),  these  higher  frequency  components  can  be 
effectively  removed  without  impairing  the  signal  level  by  lowpass  filtering.  As  with  the 
bandpass  filter,  removal  of  out  of  band  noise  will  increase  the  SNR  of  the  processed  time 
domain  signal.  The  filter  was  designed  to  pass  frequency  less  than  12  Hz,  including  the 
1 1 .4  Hz.  envelope  of  the  demodulated  signal.  The  program  allows  selection  of  the  sample 
rate  from  an  on-screen  menu.  Figure  4.5  displays  the  magnitude  and  phase  respoi  >e  of 
the  lowpass  FIR  filter  superimposed  on  the  baseband  envelope  of  the  m-sequence. 
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Figure  4.5:  Magnitude  and  Phase  Response  of  20th  Order  Lowpass  Filter. 

Figure  4.6  shows  the  same  1  second  portion  of  demodulator  output  as  Figure  4.4  but  after 
processing  through  the  lowpass  filter.  Note  that  this  figure  represents  only  one  second 
of  a  23.368  second  sequence  which  is  comprised  of  a  sequence  of  255  ones  and  zeroes. 
Additionally,  note  that  Figures  4.4  and  4.6  show  only  the  real  pan  of  the  signal.  Both 
the  real  and  imaginary  demodulated  streams  are  output  for  simultaneous  subsequent 
processing. 
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Figure  4.6:  Lowpass  Filtered,  Demodulated  255  Test  Data. 


7.  Sequence  Removal 


Both  lowpass  filtered  outputs  of  the  demodulator  are  piped  into  the  sequence 
removal  function,  finet.c.  The  data  is  passed  in  one  second  buffers  into  a  storage  array 
until  an  entire  m-sequence  has  been  accumulated  for  processing.  The  program  then 
performs  the  Fast  Hadamard  Transform  on  the  40  or  80  redundant  and  interleaved  m- 
sequences  (both  real  and  imaginary)  as  discussed  in  Chapter  HI.  The  real  and  imaginary 
components  are  then  combined  to  produce  the  magnitude  squared  (intensity)  of  the 
processed  signal.  Figure  4.7  shows  the  sequence  remover  output  for  sequence  number 
58  of  the  255  m-sequence  test  data.  A  closeup  view  of  the  same  sequence  illustrates  the 
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(b)  Closeup  of  Sequence  Remover  Output 


Figure  4.7:  Sequence  Remover  Output  for  Sequence  #58  of  255  M-Sequence  Test 
Data. 
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contribution  of  the  40  sample  points  during  the  0.087  second  correlated  digit.  The 
effectiveness  of  the  sequence  remover  is  dependent  on  the  correct  demodulation  frequency 
as  discussed  in  Chapter  m.  The  variation  of  the  audio  cassette  tape  speed  as  the  test  data 
was  digitized  created  a  Doppler-like  variation  of  signal  frequency.  Figure  4.8(a)  is  a  plot 
of  the  peak  sequence  remover  output  values  for  66  consecutive  sequences  showing  the 
variation  in  amplitude  with  time.  Figure  4.8(b)  shows  the  deviation  from  demodulation 
frequency  (57.003  Hz)  which  can  be  visually  correlated  with  the  output  magnitude  in 
Figure  4.8(a).  As  the  actual  carrier  frequency  deviates  from  the  demodulation  frequency, 
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Figure  4.8:  Frequency  Deviation  and  Sequence  Remover  Output  vs  Time. 
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the  energy  which  would  ideally  be  correlated  into  a  single  digit  is  spread  across  the  time 
series  as  noise,  which  is  illustrated  in  Figures  4.9(a  and  b).  Note  that  time  1297  seconds 
(sequence  #58  in  Figure  4.7)  had  a  deviation  of  about  0.015  Hz  (equivalent  to  .395  m/s 
relative  motion  )  while  time  962  seconds  (sequence  #43  in  Figure  4.9)  had  a  deviation  of 
over  0.03  Hz  (0.789  m/s  relative  motion). 
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Figure  4.9:  Sequence  Remover  Output  for  255  Test  Data  (Sequence  #43). 


(b)  Closcup  of  Sequence  Remover  Output 


Time  (seconds) 


70 


8.  Averaging  and  Display 

The  C  function  spewrit.c  performs  several  analytical  functions  on  the  data. 
It  scans  the  30  frequency  bins  centered  on  57.0  Hz  to  select  die  bin  with  the  highest 
amplitude.  This  corresponds  to  a  frequency  scan  from  56.903  to  57.097  Hz  which  is  the 
Doppler  frequency  shift  resulting  from  array  velocities  along  the  propagation  path  equal 
to  0  +/-  5.06  knots.  The  program  writes  an  output  file  containing  the  peak  amplitude  of 
the  scanned  bins,  the  corresponding  frequency  of  the  peak  bin,  and  the  average  noise  level 
in  a  24.8  Hz  band  centered  at  57.0  Hz.  Note  that  if  the  carrier  signal  level  is  lower  than 
a  noise  bin  in  this  band  an  incorrect  value  will  be  selected  which  will  produce  errors  in 
both  the  carrier  frequency  and  carrier  amplitude  estimation.  Therefore,  median  values  are 
used  in  our  estimates  because  the  median  is  less  sensitive  to  outliers  than  the  mean. 
These  output  files  are  in  the  form  of  a  channel-time  matrix  and  can  be  easily  used  to 
estimate  the  signal-to-noise  ratio  and  frequency  stability  of  the  carrier  wave. 

Figure  4.10  is  a  gray-scale  plot  of  the  same  66  255-digit  m-sequences  of  the 
test  data,  vertically  stacked  so  that  corresponding  digits  of  separate  sequences  are 
aligned.  When  a  signal  level  is  much  smaller  than  the  ambient  noise  level,  as  it  was  for 
the  Heard  Island  Feasibility  Experiment,  it  is  often  advantageous  to  coherently  average 
several  sequences  or  FFT’s  so  that  random  uncorrelated  noise  will  tend  to  cancel  and 
correlated  signal  energy  will  stand  out.  The  C  function,  writer. c,  writes  the  repetitive  re¬ 
sequence  outputs  of  the  sequence  removal  module  or  the  repetitive  autospectral  density 
outputs  of  the  power  spectrum  module  to  a  file  such  th3t  corresponding  tim^s  or 
frequencies  are  aligned.  The  values  in  the  output  file  of  writer.c  can  then  be  easily 
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plotted  as  a  three  dimensional  surface,  averaged,  and  plotted  as  a  two  dimensional  subset 
of  the  file. 
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Figure  4.10:  Gray-Scale  Plot  of  Sequence  Remover  Output  for  66  255-digit  M- 
Sequences  of  Test  Data. 
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C.  PRELIMINARY  DATA  ANALYSIS 


1.  Objective 

The  goals  of  the  preliminary  data  analysis  was  to  determine  the  signal-to-noise 
ratio  of  the  received  signal  and  determine  whether  it  was  sufficient  (greater  than  18  dB 
post-processing)  to  permit  arrival  time  estimation  within  the  prescribed  10  millisecond 
accuracy  for  use  in  a  long  term  global  warming  experiment.  Additionally,  sequence 
removal  was  conducted  in  an  effort  to  directly  analyze  arrival  structure  and  stability.  In 
particular,  we  desired  to  find  the  maximum  coherent  integration  time  of  the  reception 
which,  in  turn,  would  determine  the  maximum  processing  gain. 

2.  Frequency  Domain  Processing 
a.  Autospectral  Density 

Autospectral  density  (or  power  spectrum)  estimation  was  performed  using 
the  modified  parallelogram  method  [Reference  28].  The  process  yields  the  mean  square 
intensity  or  power  of  the  signal  as  a  function  of  frequency  and  is  available  as  a  packaged 
module  of  Lab  Workbench.  All  frequency  domain  processing  utilized  the  maximum  Fast 
Fourier  Transform  input  points  (N=32768),  Hamming  windowing,  and  50%  overlap  of 
adjacent  FFT’s.  The  sample  interval  (At)  was  4.386  and  2.193  milli-seconds  for  the  228 
Hz  and  456  Hz  sample  rates  respectively  and  the  corresponding  FFT  block  lengths,  given 
by  (N-l)  At,  were  143.715  seconds  and  71.857  seconds  respectively.  The  frequency 
resolution  or  binwidth  of  the  power  spectrum,  (B,)  given  by  1/A,  is  then  6.958  and  1  '.916 
millihertz,  again  dependent  on  the  sample  rate. 
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Analysis  in  the  frequency  domain  provided  a  quick  assessment  of 
equipment  operation  as  well  as  an  indication  of  relative  signal  and  noise  levels. 
Performance  of  individual  channels  can  be  estimated  from  the  deviation  of  their  output 
from  the  median  levels.  Channels  14-19,  21-23,  28,  29,  and  31  had  indications  of 
extremely  high  or  low  output.  Only  channels  whose  output  was  consistently  within  3 
standard  deviations  of  the  array  median  were  used  to  estimate  signal  and  noise  values. 
Figure  4. 1 1  shows  the  median  peak  spectral  value  (of  all  valid  channels)  for  the  01270222 


Channels  0-13,  20,  24,  27,  30  only 


Median 
1.1 803e-l  3 


1000  2000  3000  4000  5000  6000 

Time  (seconds) 

Figure  4.11:  Median  57.0  Hz  Spectral  Peak  of  01270222  CW  Signal. 

CW  signal  in  a  0.193  Hz  band  centered  at  57.0  Hz  (corresponding  to  +/-  5  kts  Doppler) 

as  a  function  of  time.  The  median  channel  spectral  noise  value  in  a  band  from  45  to  69 

Hz  is  plotted  in  Figure  4.12.  Although  band  noise  is  not  strictly  applicable  for  a  CW 

signal,  using  a  24  Hz  band  for  noise  summing  provides  a  reasonable  estimation  of 

expected  noise  levels  for  the  m-sequences.  The  noise  value  was  obtained  for  each  power 

spectrum  by  summing  the  values  in  each  frequency  bin  of  the  24  Hz  band  and  then 

dividing  by  the  number  of  bins  (862).  Finally,  the  median  value  of  all  valid  channel  -  was 
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Figure  4.12:  Median  12  Hz  Band  Noise  of  01270222  CW  vs  Time. 

plotted  as  a  function  of  time.  Taking  the  median  again  (of  the  median  channel  values) 

provides  an  overall  median  value  of  1.18xl0'13  volts2/Hz  for  the  signal  level  and  6.39x10' 

11  volts2/Hz  for  the  noise  value.  These  values  can  be  converted  to  acoustic  power  levels 

of  19  dB  rel  IpPa  for  the  signal  and  46.4  dB  rel  IpPa  for  the  noise,  as  discussed  below. 

Transmission  loss  over  the  mean  acoustic  path  would  then  be  about  213  dB  -  19  dB  = 

194  dB. 

b.  Signal-to-N oise  Ratio 

The  rms  power  in  a  given  frequency  bin  is  obtained  by  multiplying  the 
autospectral  density  G„  times  the  frequency  binwidth,  (B,).  Since  transducer  output 
voltages  were  recorded,  the  measured  rms  power  is  in  units  of  volts  squared  per  Hertz  and 
can  be  related  to  acoustic  power  by  the  sensitivity  of  the  hydrophones  (-170  Db  rel  IpPa 
or  3.162  x  10*  pPa/V).  Acoustic  power  is  the  pressure  force  times  the  velocity  of  the 
pressure  wave  (which  is  in  phase  with  the  force  at  large  distances  from  the  source). 
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Therefore  output  of  the  autospectral  density  estimator  can  be  converted  directly  to 
acoustic  power  (P)  if  desired  by 


Determining  the  SNR  of  the  CW  signals  is  simply  a  matter  of  measuring 
the  signal  power  in  the  correct  Doppler  shifted  carrier  frequency  bin  and  then  dividing 
the  result  by  the  total  noise  power  in  the  FFT  band  (summing  the  noise  power  in  all  bins). 
This  is  accomplished  by  simply  dividing  the  channel  median  of  the  spectral  peak  matrix 
(Figure  4.11)  by  the  channel  median  of  the  band  noise  matrix  (Figure  4.12).  Figure  4.13 
illustrates  the  result.  The  overall  median  signal-to-noise  ratio  for  the  reception  can  then 
be  calculated  as  10  log  (median  signal  /  median  noise)  which  equals  -27.6  dB.  Recall 
that  the  total  noise  is  significantly  reduced  by  clipping  and  filtering,  increasing  die  SNR. 


Time  (seconds) 


Figure  4.13:  Median  SNR  of  01270222  CW  Signal. 
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c.  Arrival  Structure 


Continuous  wave  signals  contain  100%  of  the  power  in  the  carrier  and 
are  therefore  the  optimum  waveform  for  resolving  the  spatial  distribution  of  energy  over 
the  water  column.  To  analyze  the  arrival  structure  of  the  signal  the  matrix  of  carrier  bin 
autospectral  density  values  was  formed.  The  matrix  dimensions  are  32  channels  by  N 
FFT  periods,  where  N  is  the  number  of  FFT  periods  in  the  reception,  so  that  each  matrix 
element  contains  the  autospectral  density  for  the  corresponding  channel  and  FFT  period. 
This  single  arrival  structure  matrix  contains  both  the  temporal  and  spatial  variability 
information  for  a  single  reception.  The  median  value  of  the  18  valid  channels,  plotted 
in  Figures  4.11  through  4.13,  provides  the  temporal  arrival  pattern  of  the  power  level 
during  the  entire  reception.  A  median  value  of  the  N  FFT’s  produces  a  plot  of  spatial 
variability  in  signal  power  level  over  the  depth  of  the  array.  The  time  averaged  FFT  is 
shown  in  Figure  4.14.  When  compared  to  the  sound  velocity  profile,  this  spatial 
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Figure  4.14:  Median  SNR  (in  24  Hz  band)  vs  Channel  for  01270222  CW 
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distribution  of  energy  might  provide  a  clue  to  the  predominant  modal  composition  of  the 
signal  and  thus  the  primary  propagation  means. 

The  median  channel  frequency  of  the  01270222  CW  signal  is  indicated 
in  Figure  4.15.  This  figure  demonstrates  that  the  frequency  was  remarkably  constant 
considering  that  both  arrays  were  subject  to  motion.  Again,  the  large  deviations  from  the 
median  are  probably  attributable  to  inaccurate  spectral  peak  selection. 
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Figure  4.15:  Median  Frequency  of  01270222  CW  Signal. 

3.  Time  Domain  Processing 

A  pessimistic  estimate  of  the  time  domain  signal-to-noise  ratio  is  made  by 
finding  the  mean  amplitude  of  all  the  points  in  a  particular  sequence,  not  trying  to  sort 
out  signal  from  noise.  The  peak  magnitude  can  then  be  divided  by  this  value  to  obtain  an 
estimate  of  the  signal-to-noise  ratio.  However,  if  the  post-processed  signal  has 
insufficient  amplitude  to  be  reliably  detected  in  the  noise  field,  the  effectiveness  of  peak- 
picki  ig  is  diminished.  Unfortunately,  this  was  the  case  for  most  m-sequences  analyzed 
in  this  thesis.  Although  there  were  peaks  in  the  data  wliicli  were  stationary  in  time 
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between  several  sequences,  they  rarely  had  sufficient  amplitude  or  duration  to  distinguish 
them  from  similar  features  in  the  noise  field.  In  addition,  if  multiple  arrivals  are  present 
within  the  same  sequence,  the  difficulty  in  discriminating  the  signal  from  noise  will  be 
increased.  Note  also  that  cross-channel  average  values  cannot  be  used  effectively  in  the 
time  domain  because  of  possible  phase  differences  between  hydrophones. 

The  2047-digit  m-sequence  which  was  received  at  01270005  GMT  provided 
maximum  gain  and  a  relatively  stable  frequency.  Figure  4.16(a  and  b)  depicts  the  highest 


peak  detected  in  this  reception  that  could  be  reasonably  associated  with  the  signal.  The 
peak  was  9  dB  above  the  median  noise  level  in  the  180  second  sequence  and  occurred 
at  112.6  seconds  from  the  start  of  the  7th  sequence.  Figure  4.16(b)  is  a  close-up  view  of 
the  same  peak,  depicting  the  width  of  one  m-sequence  digit  (.0877  seconds  or  20  sample 
points).  This  particular  peak,  although  higher  than  any  others,  was  not  repeated  in  other 
sequences  at  a  corresponding  time  from  the  sequence  start.  Occurring  only  2  hours  before 
the  01270222CW  signal,  the  01270005  reception  had  similar  noise  levels.  The  coded 
signal  level  is  3  dB  less  than  the  CW  level  and  thus  the  expected  pre-processed  SNR  is 
approximately  -30  dB.  Since  the  maximum  gain  of  the  2047  digit  pulse  compression  is 
33  dB,  we  can  only  expect  a  maximum  output  SNR  of  about  3  dB  under  ideal  conditions. 

If  the  maximum  coherent  integration  time  of  the  arrival  exceeds  the  length  of 
two  consecutive  sequences  (359  seconds  in  this  case),  additional  gain  can  be  achieved  by 
coherently  summing  adjacent  sequences.  In  addition,  more  credence  can  be  given  to  an 
arrival  peak  that  is  repetitive.  Sequence  averaging  was  performed  in  3  minute  increments 
from  6  minutes  up  to  one  hour  in  an  effort  to  locate  stable  peaks.  A  peak  was  found  in 
sequence  number  8  of  channel  5  at  169.25  seconds  from  the  beginning  of  the  sequence 
which  also  occurred  in  other  sequences.  Figure  4.17(a)  and  4.17(b)  illustrate  this  peak. 
Note  that  this  peak  not  only  has  a  smaller  amplitude  than  Figure  4.16  but  is  also 
narrower,  similar  to  the  test  data  case  in  Figure  4.8.  The  peak  was  detectable  above  the 
noise  in  sequences  7,  8,  9,  10,  and  13.  Figure  4.18  shows  the  mean  amplitude  of  this 
peak  averaged  over  sequences  6  through  13.  This  mean  peak  is  thus  a  coherent  avc  rage 
between  8  sequences  (about  24  minutes)  and  is  roughly  6.7  dB  above  the  mean  noise. 
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(a)  Sequence  #8 


Time  Since  Start  of  Sequence  (sec) 


Time  Since  Start  of  Sequence  (sec) 


Figure  4.17:  A  Persistent  Peak  in  Sequence  8  of  01270005  2047-Digit  M-Sequence. 


Time  Since  Start  of  Sequence 


Figure  4.18:  Mean  Amplitude  of  a  Persistent  Peak  in  01270005  2047-Digit  M- 
Sequence. 
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Further  evidence  that  the  peak  in  Figures  4.17  and  4.18  is  actually  an  arrival 
of  2047  m-sequence  code  is  that  they  were  not  detected  at  all  in  sequences  1  through  6. 
The  predicted  first  arrival  of  each  sequence  was  3  hours  and  20  minutes  after  the 
transmission  commenced  due  to  acoustic  propagation  at  a  nominal  1500  m/s  over  18  Mm. 
Since  this  particular  recording  commenced  3  hours  and  5  minutes  after  the  transmission, 
the  first  peak  would  be  expected  to  occur  15  minutes  into  the  recording  after  sequence 
number  5.  Figure  4.19  is  a  plot  of  the  mean  level  in  a  0.0877  second  band  around  the 
peak  as  a  function  of  sequence  number.  The  energy  from  the  peak  can  be  seen  to 
increase  after  sequence  number  5,  as  expected. 
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Figure  4.19:  Mean  Amplitude  of  01270005  Peak  vs  Sequence  Number. 

The  stability  of  an  arrival  determines  its  maximum  coherent  integration  time 
since  any  time  shift  due  to  Doppler  effects  or  dispersion,  or  amplitude  decrease  will 
diminish  summing  gains.  A  three  dimensional  perspective  of  the  stability  of  the  peak  in 
Figure  4.17  is  shown  in  Figure  4.20.  Note  that  the  peak  varies  in  both  amplitude  as  a 
function  of  sequence  number  and  time  position  within  any  given  sequence.  In  addition 
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Sequence  Number 


to  a  myriad  of  oceanic  causes,  the  variation  could  be  produced  or  compounded  by 
deviation  of  the  received  carrier  frequency  from  the  demodulation  frequency,  as  discussed 
in  Section  4.B. 7.  A  grayscale  plot  of  this  same  reception  is  shown  in  Figure 
4.21. 


SEQUENCE  REMOVER  OUTPUT  01270005 


Time  Since  Start  of  Sequence 


Figure  4.20:  Variation  of  Amplitude  and  Position  of  2047-Digit  M-Sequence. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  SUMMARY  OF  RESULTS 

The  Naval  Postgraduate  School  in  conjunction  with  The  Monterey  Bay  Aquarium 
Institute,  Woods  Hole  Oceanographic  Institution,  and  Massachusetts  Institute  of 
Technology  successfully  participated  in  1991  The  Heard  Island  Feasibility  Experiment, 
demonstrating  the  ability  to  receive  global  transmission  of  acoustic  signals.  The  18 
Megameter  transmission  may  have  been  the  furthermost  reception  ever  of  an 
electronically  generated  acoustic  signal  in  an  oceanic  medium.  Over  1.5  gigabytes  of  data 
were  collected  from  a  32  hydrophone  array  which  intersected  the  deep  sound  channel  in 
a  horizontal  zone  of  optimum  reception. 

The  Naval  Postgraduate  School  has  gained  a  real-time  processing  capability  for 
tomographic  m-sequence  signals  using  the  Laboratory  Workbench  environment  of  the 
Concurrent  6400  data  acquisition  system.  The  potential  to  process  phase-encoded  signals 
as  they  are  received  will  enhance  the  schools  efforts  in  future  acoustic  tomography  efforts. 
Adjustments  to  experimental  parameters  can  now  be  made  in  a  timely  manner  to  optimize 
the  result. 

Only  a  very  preliminary  analysis  of  the  Heard  Island  data  has  been  performed  in 
order  to  provide  an  overview  of  signal  levels,  integration  times,  and  stability.  Detection 
of  the  Heard  Island  was  possible  in  both  the  frequency  domain  and  the  time  domain. 
Spectral  analysis  of  CW  signals  indicated  a  relatively  stable  carrier  frequency  (+/-  about 
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0.01  Hz  corresponding  to  about  0.5  kts  of  relative  motion).  The  CW  signal  had  a  pre- 
processed  level  of  about  -28  +/-  3  dB  which  would  imply  an  approximate  -31  dB  level 
for  a  coded  signal  in  a  similar  noise  field. 

Time  domain  processing  using  pulse  compression  of  the  coded  m-sequence  enabled 
identification  and  resolution  of  occasional  arrivals  of  the  2047-digit  m-sequence. 
However  the  post-processed  signal  level  was  only  about  9  dB  for  the  maximum  peak 
found  and  most  arrival  peaks  had  characteristics  so  similar  to  the  noise  field  that  reliable 
identification  was  not  possible.  The  maximum  coherent  integration  time  found  was  7 
sequences  or  about  21  minutes  of  the  2047-digit  sequence.  The  pre-processed  levels  were 
so  low  that  small  increases  in  the  ambient  noise  level  cause  the  post-processed  SNR  to 
fall  below  zero  and  prevent  detection.  Detection  of  the  signal  in  both  the  frequency  and 
time  domains  was  intermittent  as  the  signal  often  fell  below  the  detection  threshold. 

B.  CONCLUSIONS 

I.  Signal  Processing 

The  signal  processing  software  functions  satisfactorily  as  designed  but  is 
sensitive  to  deviations  of  the  carrier  from  the  demodulation  frequency,  which  occur  as  a 
result  of  relative  motion  between  the  source  and  receiver.  A  deviation  of  .01  Hz  lowered 
the  output  amplitude  by  about  1.5  dB.  A  frequency  deviation  of  more  than  0.02  Hz  is 
generally  sufficient  to  make  the  signal  undetectable  in  the  time  domain.  Resampling  of 
the  signal  may  be  necessary  to  minimize  the  effects  of  large  Doppler  frequency  shifts 
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during  moving  ship  tomography  experiments  with  low  SNR.  The  signal  processing 
objectives  of  this  thesis,  as  outlined  in  Chapter  I  have  been  fully  met. 

2.  Heard  Island  Feasibility  Experiment 

Analysis  of  the  Heard  Island  data  has  not  been  completed.  The  very  cursory 
processing  conducted  in  thesis  was  performed  in  order  to  verify  the  validity  of  the 
processing  procedures  and  to  begin  to  quantify  gross  characteristics  of  the  data.  Both 
frequency  and  time  domain  processing  are  viable  analysis  tools  for  this  data,  each 
possessing  its  own  advantages. 

The  SNR  calculations  in  Chapter  IV  using  the  median  signal  and  noise  levels 
indicate  that  the  maximum  processing  gain  of  the  longest  sequence  used  (33  dB)  is  barely 
sufficient  to  make  the  signal  detectable  in  the  noise  field.  The  average  pre-processed 
SNR  of  -28  dB  is  indicative  of  a  19  dB  acoustic  power  level  rel  IpPa  or  in  other  words 
a  194  dB  transmission  loss.  Indeed,  the  processed  signal  is  tenuous,  at  best,  and  often 
undetectable  in  the  noise  field.  Additional  gain  is  required  in  order  to  reliably  identify 
and  analyze  the  arrival  structure. 

The  maximum  observed  processed  SNR  in  the  time  domain  was  9  dB  which 
corresponds  to  an  arrival  uncertainty  of  29  milliseconds,  using  Eq.  3.1.  This  is 
approximately  3  times  the  maximum  uncertainty  objective  of  the  Heard  Island  Feasibility 
Experiment.  This  maximum  peak  was  not  repetitive  and  was  thus  not  reliably  identifiable 
or  conducive  to  averaging  gain.  Arrival  peaks  had  pulse  widths  on  the  order  of  0.1 
seconds.  Typical  peaks  had  levels  on  the  order  of  4  dB  which  corresponds  to  an  arrival 
time  uncertainty  of  about  44  milliseconds.  Therefore,  without  further  processing  gai'i,  the 
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propagation  from  Heard  Island  to  Monterey  exhibits  too  much  transmission  loss  to  meet 


the  experiment  objectives. 


C.  RECOMMENDATIONS  FOR  ADDITIONAL  WORK 

There  is  much  work  left  to  do  in  analysis  of  the  Heard  Island  data  and  improving 
the  acoustic  tomography  capability  at  Naval  Postgraduate  School.  Some  recommendations 
in  this  regard  are  listed  below: 

1 .  Additional  software  development  should  be  pursued  to  facilitate  real-time  Doppler 
frequency  search  and  compensation.  Automated  phase  correction  would  maximize  the 
efficiency  of  the  entire  sequence  removal  process.  The  effects  of  relative  motion  are 
significant  and  must  be  removed  in  order  to  have  persistent,  high  amplitude  arrival 
peaks. 

2.  Continued  efforts  in  reducing  the  noise  in  the  recorded  data  would  improve  the 
output  SNR.  Pre-whitening  of  the  noise  using  standard  methods  such  as  an  auto¬ 
regressive  filter  may  improve  the  performance  of  the  noise  filters  as  well  as  enhance 
coherent  averaging.  Also,  the  origin  of  the  large  noise  transients  that  occurred  in  all 
channels  should  be  investigated. 

3.  Array  beamforming  can  be  implemented  to  obtain  up  to  15  dB  gain  ideally  for  all 
32  hydrophones.  Up  to  12  dB  of  gain  could  be  achieved  for  the  18  functional 
hydrophones.  In  addition,  modal  beamforming  would  act  as  a  modal  filter  and  allow 
discrimination  of  the  predominant  modes. 

4.  Complete  the  initial  survey  of  the  data  by  creating  matrices  of  the  output  data  and 
identifying  signal  peaks.  If  sufficient  gain  can  be  realized,  a  detailed  study  of  the 
arrival  structure  should  be  undertaken  to  begin  to  quantify  the  effects  of  dispersion  and 
oceanic  influences  on  the  various  propagation  paths.  Within  the  acoustic  data  from 
this  experiment  lies  the  signature  of  every  oceanic  process  between  Heard  Island  and 
Monterey  for  a  period  of  7  days.  Further  analysis  could  provide  information  on  the 
variability  of  the  marine  environment  over  this  period.  Again,  Doppler  effects  must 
be  removed  since  they  would  mask  any  short  term  oceanographic  effects  on 
propagation. 

5.  The  Heard  Island  Feasibility  Experiment  did  not  utilize  precise  timing  between  the 
source  and  the  receivers  because  the  objective  was  determining  the  accuracy  of  the 
arrival  time  and  not  the  arrival  time  itself.  For  future  tomography  work,  a  precision 
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time  synchronization  system  should  be  devised  which  could  supplement  the  acoustic 
data  with  clocking  information. 

6.  The  volume  of  stored  data  in  a  tomography  experiment  can  be  enormous. 
However,  since  each  digit  is  sampled  numerous  times,  there  is  a  redundancy  in  the 
data  that  can  be  reduced  by  decimation  to  conserve  storage  space.  The  C  function 
fmetl.c  has  the  ground  work  for  a  decimation  scheme  already  installed  in  the  form  of 
a  decimation  factor  that  is  menu-selectable.  This  feature  should  be  expanded  upon  in 
future  experiments. 

7.  Due  to  the  oversampling  described  above,  the  top  of  the  correlation  peaks  in  the 
m-sequence  are  flat  vice  ’pointy’  as  would  be  the  case  for  a  single  sample  per  digit. 
To  more  accurately  resolve  the  actual  arrival  time,  it  would  be  advantageous  to 
implement  a  method  of  locating  a  particular  feature  of  die  peak  such  as  its  center  or 
forward  edge.  The  center  can  be  resolved  by  cross-correlating  the  peak  with  a  square 
pulse  of  the  same  width  which  results  in  a  ’pointed’  peak. 


Finally,  since  the  propagation  path  from  Heard  Island  to  Monterey  exhibits  too 
much  transmission  loss  to  be  useful  for  a  long  term  study  of  global  warming,  it  would 
be  useful  to  conduct  ray  tracing  between  alternate  sites  in  the  Pacific  basin.  The 
Hamiltonian  methods  discussed  in  Chapter  HI  could  be  used  to  determine  the  optimum 
locations  of  both  a  source  and  receiver  to  perform  long  term  monitoring  of  global 
temperatures. 
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APPENDIX  A 


SCHEDULE  OF  RECORDED  DATA 


GREENWICH 

MEAN 

TIME 

LOCAL 

TIME 

(PST) 

SIGNAL 

TYPE 

DATA 

LENGTH 

(MIN) 

DATA 

SIZE 

(MBytes) 

SAMPLE 

FREQ. 

(Hz) 

01261525 

01260725 

CW 

121 

106 

228  Hz 

01261800 

01261000 

3/10 

107 

93 

228  Hz 

01262058 

01261258 

255  M 

38 

33 

228  Hz 

01270005 

01261605 

2047  M 

73 

64 

228  Hz 

01270222 

01261922 

CW 

85 

75 

228  Hz 

10270608 

01262208 

3/10 

75 

66 

228  Hz 

01271207 

01270407 

511  M 

80 

70 

228  Hz 

01271505 

01270705 

CW 

158 

138 

228  Hz 

01271807 

01271007 

3/10 

153 

134 

228  Hz 

01292115 

01291315 

255  M 

138 

241 

456  Hz 

01300011 

01291611 

1023  M 

72 

126 

456  Hz 

01312108 

01311308 

255  M 

28 

49 

456  Hz 

02010011 

01311611 

2047  M 

18 

31 

456  Hz 

Notes:  1.  GMT(recdved)  =  GMT(IX)  +  3hrs  &  20m in 

2.  PST(received)  =  GMT(W  -  8  hrs  +  3hrs  &  20m in 

3.  Acoustic  Travel  Time  From  Heard  Island  to  Monterey  is  about  3  hours 

and  20  minutes. 
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Appendix  B 

SIGNAL  PROCESSING  PROGRAMS 
The  following  programs  are  listed  in  this  appendix: 

A.  RESAMP.C  -  Resampling  for  Doppler  Correction  and  Clipping. 

B.  BANDPASS  C  -  20th  order  bandpass  filter  (45-69  Hz). 

C.  FDEMOD1.C  -  Demodulator. 

D.  LOP  ASS  .C  -  20,h  order  lowpass  filter  (0-12  Hz). 

E.  FMET1.C  -  Sequence  Removal. 

F.  SPEWRIT.C  -  Writes  spectral  peak,  frequency,  and  band  noise  files. 
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RES AMP. C 

This  program  forms  a  Lab  Workbench  module  which  resamples  an  input  data 
stream  at  exactly  four  times  a  menu  specified  carrier  frequency.  It  was 
designed  to  be  used  with  a  buffer  size  of  10  seconds  so  the  header  file  of 
the  data  will  have  to  be  adjusted  accordingly.  The  program  also 
calculates  the  mean  and  standard  deviation  of  each  input  buffer  and  will 
clip  all  data  in  the  buffer  to  a  menu  specified  factor  times  the  standard 
deviation  from  the  mean  (default  is  2). 

************************************************************************* f 

static  char  SccsID[]  =  "@(#)uftest.c8. 1  (MASSOOMP)  2/5/90"; 

#include  "Iwbufunc.h" 

#include  "errors. h" 

#include  "math.h" 
float  ’idata,  ’odata; 

double  nfc=57.0,  clip=2.0,  pres [2280] [32] ,  past [2280] [32] .  futr[2280] [32] ; 
int  bufcount,  omsg=0, links=0,  state=0, itype=0,  otype=LWBFLOAT, 
dimens ion=0,  sizes [3] ={0,0,0} ,  count=0,  samples,  rate=228.0,  dx=0; 
/****************************************  ******************************** 
*  >>  Entry  point  from  LWB  —  must  be  first  function  defined  in  file  <<  * 

********************************************************** ************** i 


uftest(  action,  inA,  inB,  inC, 
int  ’action,  *inA,  »inB,  *inC, 

{ 

int  status; 
dx  =  ’debug; 
switch  (  ’action  ) 

{ 

case  SETUPCMD: status 
case  RUNCMD:  status 
case  STOPCMD: status  = 
case  MENUCMD: status  = 
case  DROPCMD: status  = 
default:  break; 

} 

return (  status  ) ; 


inD,  outX,  outY,  debug  ) 
*inD,  outX,  outY,  ’debug; 


=  ufsetup(  *inA,  outX  );  break; 
=  ufrun(  *inA,  outX  );  break; 

=  ufstop(  *inA,  outX  );  break; 

=  ufmenu(  *inA  );  break; 
ufdropQ ;  break; 


} 


SETUP 


/*  ******** ***************** 

ufsetup(  in,  out  ) 
int  in,  out; 

{ 

int  i; 

bufcount  =  0; 

/*  set  up  menu  label  and  initial  value 


************************ 


/ 
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if  (  ! lwblabelmenu(  1,  "Clip  Factor"  )  )  retum(  0  ); 
if  (  ! lwbupdatemenu(  1,  "%f",  clip)  )  retum(  0  ); 
if  (  ! lwblabelmenu(  2,  "Received  Freq."  )  )  retum(  0  ); 
if  (  ! lwbupdatemenu(  2,  "%f",  nfc)  )  retum(  0  ); 

/*  get  dimensionality  and  size  of  input  message  data  */ 
lwbgetparam(  in,  0,  LWBMSGDIM,  &dimension  ); 


if  (dx)  printf(  "resamp  dimension=%d\n" ,  dimension  ); 
lwbgetparam(  in,  0,  LWBMSGSIZES,  sizes  ); 

if  (dx)  printf  ("resamp  sizes=%d,%d,%d\n",  sizes  [0]  .sizes [1]  .sizes [2]) 
if  (  dimension  =  2  ) 

{ 

count  =  sizes [0]  *  sizes [1]; 
samples  =  sizes [0]; 

} 

else  if  (  dimension  =  1  ) 

{ 

count  =  sizes [0]; 
samples=  sizes[0]; 

} 

if  (  dx  )  printf(  "resamp  samples=%d\n" ,  samples); 

/*  check  input  data  type  */ 

lwbgetparam(  in,  0,  LWBMSGTYPE,  Aitype  ); 

if  (  dx  )  printf(  "resamp  itype=%d\n”,  itype  ); 

if  (  itype  !=  5  ) 

{ 

lwberror(  ERROR,  "Input  message  data  type  must  be  float"  ); 
return (  0  ) ; 

} 

/*  create  a  new  message  for  output  data  */ 
if  (  !lwbnewmsg(  Aomsg  )  )  retum(  0  ); 

/*  initialize  message  */ 

if  (  !lwbcopymsg(  in,  omsg  )  )  retum(  0  ); 
if  (  !lwbsamerange(  in,  omsg  )  )  retum(  0  ); 
if  (  !  lwbsamedomain(  in,  omsg  )  )  retum(  0  ); 

/*  create  a  data  buffer  for  output  message  */ 
if  (  !lwbnewbuf(  omsg  )  )  retum(  0  ); 
if  (  ! lwbgetdatbuf  (  omsg,  Aodata  )  )  retum(  0  ); 

/*  send  message  on  its  way  and  release  the  input  */ 
lwbwrite(  omsg,  out  ); 
lwbrelease(  in  ); 
return (  1  ) ; 

) 

!% **************************  run  it******************************/ 

ufrun(  in,  out  ) 
int  in,  out; 

{ 

register  int  i,j,k; 
int  mi ,oi ,ci ,f i ,pi ; 

double  rise,  run,  lower,  slope,  delt,otime,ntime,diff ,  mean[32], 
var[32],  sigma[32],  sum[32],  num[32],  max[32],  min[32]; 

/*  sanity  check  */ 

if  (  !in  i!  !omsg  )  return(  0  ); 

if  (  ! lwbgetparam(  omsg,  0,  LWBMSGLINKS,  Alinks  )  )  return(  0  ); 
if  (  links  >  0  ) 

{ 

lwbunget(  in  ); 
return (  1  ); 
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if  (  ! lwbgetdatbuf (  in,  Aidata  )  )  return (  0  ); 

/***  IF  THIS  IS  THE  FIRST  BUFFER  ***/ 

if(bufcount  =  0)  /*  Note:  1st  10  sec  buffer  output  is  bogus  */ 

{ 

/***  GET  THE  MEAN  AND  STD  DEV  FOR  EACH  CHANNEL  OF  THIS  BUFFER  ***/ 
for  (  j=0;  j<32;  j++) 

{ 

sumfj]  =  0; 

for(i=0;  i<samples;  i++) 

{ 

k=j+(32*i);  /*  index  of  idatati] [j]  */ 

sumfj]  =  sumfj]  +  (double) i data [k] ; 

} 

} 

for( j=0;  j <32 ;  j++) 

{ 

meanfj]  =  sumfj]  /  (double) samples; 
diff  =  0.0; 
numfj]  =  0.0; 
for(i=0;  i<samples;  i++) 

{ 

k=j+(32*i);  /*  index  of  idatafi] [j]  */ 

diff  =  (double) idatafk]  -  meanfj]; 
numfj]  =  numfj]  +  (diff  *  diff); 

} 

varfj]  =  numfj]  /  (samples-1); 
sigmafj]  =  sqrt(varfjj) ; 
maxfj]  =  meanfj]  +  (clip  *  sigmafj]); 
minfj]  =  meanfj]  -  (clip  *  sigmafj]); 


/***  READ  INPUT  INTO  PRESENT  ARRAY  &  CLIP  ***/ 
for  (  i=0;  i<samples;  i++) 

{ 

for(j=0;  j<32;  j++) 

{ 

k=j+(32*i);  /*  index  of  idatafi] [j]  */ 

if ((double) idatafk]  <  maxfj]  &&  (double) idata  >  minfj]) 

{ 

presfi][j]  =  (double) idatafk] ; 
odatafk]  =idata[k]; 

> 

else  /*  clip  if  input  >  clip  *  sigma  */ 

{ 

if (idatafk]  >=  0.0) 

{ 

pres [i]  [ j]  =  maxfj]; 
odatafk]  =  (f loat)max[ j] ; 

} 

if (idatafk]  <  0.0) 


94 


{ 

pres  [i ]  [ j]  =  min[ j] ; 
odata[k]=(f loat)min[  j] ; 

} 

} 

} 

> 

} 

/***  FOR  ALL  REMAINING  BUFFERS  ***/ 
else 
{ 

/***  READ  INPUT  INTO  FUTURE  ARRAY  ***/ 
for(i=0;  i<samples;  i++) 

{ 

for(j=0;  j<32;  j++) 

{ 

k=j+(32*i) ;  /*  index  of  idatafi] [j]  */ 

futr[i]  [j]  =  (double) i data [k] ; 

} 

} 

/***  RESAMPLE  THE  'PRESENT'  ARRAY  ***/ 
for  (  i=0;  i<samples;  i++) 

{ 

mi=(bufcount*samples)+i ;  /*  running  index  for  new  time  scale  */ 
ntime=(double)mi/(4.0*nfc) ;  /*  new  actual  time  for  current  i  */ 
oi  =  ntime  *  rate;  /*old  index  just  lower  than  newtime*/ 

ci=oi-(bufcount*samples) ;  /*  current  index  lower  than  newtime  */ 
otime=(double)oi/rate;  /*  old  actual  time  lower  than  newtime  */ 
run  =  1.0  /  rate;  /*  time  difference  between  old  samples  */ 
delt  =  ntime-otime;  /*  time  shift  for  this  sample  V 

for(j=0;  j <32 ;  j++) 

{ 

k=j+(32*i);  /*  index  of  odata[i] [j]  */ 

if  (ci  >=  0  &&  ci+1  <=  samples)  /*  both  values  in  present  */ 
{ 

lower  =  pres [ci] [ j] ; 

rise  =  presfci+l] [j]  -  lower; 

} 

else  if(ci+l  ==  0)  /*  values  split  pres/past  */ 

{ 

lower  =  past [samples-1] [ j] ; 
rise  =  pres[0][j]  -  lower; 

} 

else  if  (ci  <  0  M  ci+1  ==  0)  /*  both  values  in  past  */ 

{ 

pi  =  samples*  ci; 
lower  =  past [pi] [ j] ; 
rise  =  past [pi+1] [ j]  -  lower; 

} 

else  if  (ci  =  count-1)  /*  values  split  pres/futr  */ 
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{ 

lower  =  pres[samples-l] t j] ; 
rise  =  futr[0][j]  -  lower; 

> 

else  if  (ci  >  samples-l)  /*  both  values  in  futr  */ 

{ 

fi  =  ci  -  samples; 
lower  =  futr[f i] [j] ; 
rise  =  futrff i+1] [j]  -  lower; 

> 

slope  =  rise  /  run; 

odata[k]  =  (float) (lower  +  (slope  *  delt)); 

} 

} 

/***  GET  THE  MEAN  AND  STD  DEV  FOR  EACH  CHANNEL  OF  THE  OUTPUT  ***/ 
for  (  j=0;  j<32;  j++) 

{ 

sumfj]  =  0; 

for(i=0;  i<samples;  i++) 
k=j+(32*i);  /*  index  of  odata[i][j]  */ 

sumfj]  =  sum[j]  +  (double)odata [k] ; 

} 

for( j=0;  j  <32 ;  j++) 

{ 

meanfj]  =  sumfj]  /  (double) samples; 
diff  =  0.0; 
num[j]  =0.0; 
for(i=0;  i<samples;  i++) 

{ 

k=j+(32*i);  /*  index  of  odata[i][j]  */ 

diff  =  (double) odata[k]  -  meanfjj; 
num[j]  =  num[j]  +  (diff  *  diff); 

> 

ave[j]  =  num[j]  /  (samples-l); 
sigma [j]  =  sqrt (var [ j] ) ; 
max[j]  =  meanfj]  +  (clip  *  sigma[j]); 
min[j)  =  mean[j]  -  (clip  *  sigma[j]); 

> 

/*  CLIP  OUTPUT  */ 
for(i=0;  i<samples;  i++) 

{ 

for( j=0;  j<32;  j++) 

{ 

k=j+(32*i);  /*  index  of  odata[i][j]  */ 

if ((double) odata[k]  >  max[j]  !!  (double) odata[k]  <  min[j]) 

{ 

if (odatafk]  >=  0.0) 
odata[k]  =  (f loat)max[ j] ; 
if (odatafk]  <  0.0) 
odatafk]  =  (f loat)min{ j] ; 

} 
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} 


} 

/***  shift  the  data  arrays  toward  the  past  ***/ 
for  (  i=0;  i<samples;  i++) 

{ 

for(j=0;  j<32 ;  j++) 

{ 

past  [i]  [j]  =  pres [i]  [j] ; 
pres [i]  [j]  =  futr [i]  [j] ; 

> 

} 

} 

bufcount++; 

lwbwrite(  omsg,  out  ); 
lwbrelease(  in  ); 
return (  1  ) ; 

} 

/ft***********************  sxop  ****************»*************/ 

ufstop(  in,  out  ) 
int  in,  out; 

{ 

if  (  Jin  !!  lout  )  retum(  0  ); 
lwbwrite(  omsg,  out  ); 
lwbrelease(  in  ); 
return (  1  )  ; 

} 

/***********************  MF.NI  1  ********************************/ 

ufmenu(  key  ) 
int  key; 

{ 

int  i ; 

if  (  key  =  1  ) 

{ 

lwbgetdouble(  "Enter  Clip  Factor",  0.0,  500.0,  Aclip); 
if  (  !lwbupdatemenu(  1,  "%f",  clip)  )  retum(  0  ); 

} 

if  (  key  =  2  ) 

{ 

lwbgetdouble(  "Enter  Received  Frequency”,  0.0,  500.0,  Anfc); 
if  (  ! lwbupdatemenu(  2,  "%f",  nfc)  )  retum(  0  ); 

} 

return (  1  ); 

} 

/***********************  DROP  **********************************/ 

ufdropO 

{ 

if  (  !lwbfreemsg(  Aomsg  )  )  retum(  0  ); 
return (  1  ) ; 

} 
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/*************************************»**********#**************»******* 

BANDPASS. c 

This  function  implements  a  20th  order  Finite  Impulse  Response  (FIR) 
bandpass  filter  where  the  filter  coefficients  were  computed  using 
the  FIR1  routine  in  MATLAB.  If  W(n)  denotes  a  window,  where  l<n<N, 
and  the  impulse  response  of  the  ideal  filter  is  h(n) ,  where  h(n)  is 
the  inverse  Fourier  transform  of  the  ideal  frequency  response,  then 
the  windowed  digital  filter  coefficients  are  given  by 
b(n)=  W(n)h(n) ,  l<n<N.  The  user  can  select  either  228  or  456  hz  sample 
rate  as  appropriate  for  the  Heard  Island  Data. 

************************************************************************  j 

static  char  SccsID[]  =  "@(#)uftest.c8. 1  (MASSCOMP)  2/5/90”; 

#include  "lwbufunc.h" 

♦include  "errors. h" 
float  ‘idata,  ‘odata; 

int  omsgS,  linksS,  states,  itype=0,  otype=LWBFLOAT ,  dimensions, 
sizes  [3] ={0,0,0}  ,  counts,  rate=228,  dxS; 
double  b_228[21]  =  {9.679120351 15705eS4,  7. 59486633 167777e-l 9, 

7.34078321 026075e-03 ,  -8 . 928452 1 000 1 374e- 1 8 , 

-4.4637751 2008602e-02 ,  2.422243186181 52e- 1 7 , 

1 .21524732772206e-01 ,  -3.05358871759598e-17, 

-2 . 059 1 9 1 47 1 9 1 579e-0 1 ,  1 . 428569 1411 9550e- 1 7 , 

2 . 430909953204 1 8e-0 1 ,  1 . 428569 1411 9550e- 1 7 , 

-2.05919147191 579e-0 1,-3. 0535887 1 759598e- 1 7 , 

1.21 524732772206e~0 1 ,  2 . 422243 186181 52e- 1 7 . 

-4 . 463775 1 2008602e-02 ,  -8 . 928452 1 000 1 374e- 1 8 , 

7 . 3407832 1 026076e-03 ,  7 . 59486633 1 67777e-l 9 , 

9.679120351 15705e-04> , 

b_456  [21]  =  {3. 01 3898085 17666eS7,  9. 658849994529 10e-03, 

2 . 44744230763637e-02 ,  3.0011 45454787 1 5e-02 , 

-1 .259082751 56001 e-06,  -6. 760841 2358034 le-02, 

-1 . 26039643870063e-01 ,  -1 . 0938997346 1736e-01 , 

1.11 960800803549e-06 ,  1 . 36893 1 3 1 470553e-0 1 , 

1 .98959837606674e-01 ,  1 .368931 3 1470553e-01 , 

1 . 1 1 960800803549e-06 ,  -1 .0938997346 1736e-01 , 

-1 . 26039643870063e-0 1 ,  -6. 760841 2358034 le-02, 

- 1 . 25908275 1 5600 1 e-06 ,  3.0011 45454787 1 5e-02 , 

2 . 447442307 63637e-02 ,  9 . 6588499945291 0eS3 , 

3.01 3898085 1 7666e-07} , 

stor [20]={1 ,0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0,0}, 

B [2 1 ] ,  otemp[500]; 

/************************************************************************ 

*  »  Entry  point  from  LWB  —  must  be  first  function  defined  in  file  <<  * 
******************************** **************************************** / 

uftest(  action,  inA,  inB,  inC,  inD,  outX,  outY,  debug  ) 
int  ‘action,  *inA,  *inB,  *inC,  *inD,  outX,  outY,  ‘debug; 

{ 

int  status; 
dx  =  ‘debug; 
switch  (  ‘action  ) 

{ 
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case  SETUPCMD: status  =  ufsetup(  *inA,  outX  );  break; 
case  RUNCMD:  status  =  ufrun(  *inA,  outX  );  break; 
case  STOPCMD: status  =  ufstop(  *inA,  outX  );  break; 
case  MENUCMD: status  =  ufrnenu!  *inA  );  break; 
case  DROPCMD: status  =  ufdropO ;  break; 
default:  break; 

} 

return (  status  ) ; 

} 

ufsetup(  in,  out  ) 
int  in,  out; 

{ 

int  i ; 

/*  set  up  menu  label  and  initial  value  */ 

if  (  ! lwblabelmenu(  1,  "Sample  Rate"  )  )  retum(  0  ); 
if  (  !  lwbupdatemenu(  1,  "%d",  rate)  )  retum(  0  ); 
if  (rate=228) 

for(i=0;  i  <  21;  i++)  B[i]  =  b_228[i]: 
if  (rate=456) 

for(i=0;  i  <  21;  i++)  B[i]  =b_4S6[i]; 

/*  get  dimensionality  and  size  of  input  message  data  */ 
lwbgetparam(  in,  0,  LWBMSGDIM,  Adimension  ); 
if  (  dx  )  printf(  "dimension=%d\n" ,  dimension  ); 
lwbgetparam(  in,  0,  LWBMSGSIZES,  sizes  ); 

if  (  dx  )  printf(  "sizes=%d,%d,%d\n" ,  sizes [0],  sizes [1],  sizes [2]  ) 
if  (  dimension  —  2  ) 

count  =  sizes [0]  *  sizes[l]; 
else  if  (  dimension  =  1  ) 
count  =  sizes [0] ; 

if  (  dx  )  printf(  "count=%d\n" ,  count  ); 

/*  check  input  data  type  */ 

lwbgetparam(  in,  0,  LWBMSGTYPE,  Aitype  ); 
if  (  dx  )  printf(  " i type=%d\n" ,  itype  ); 
if  (  itype  !=  5  ) 

{ 

lwberror(  ERROR,  "Input  message  data  type  must  be  float"  ); 
return (  0  ) ; 

} 

/*  create  a  new  message  for  output  data  */ 
if  (  !lwbnewmsg(  Aomsg  )  )  return (  0  ); 

/*  initialize  message  */ 

if  (  !lwbcopymsg(  in,  omsg  )  )  retum(  0  ); 
if  (  ! lwbsamerange(  in,  omsg  )  )  retum(  0  ); 
if  (  ! lwbsamedomain(  in,  omsg  )  )  return!  0  ); 

/*  create  a  data  buffer  for  output  message  */ 
if  (  !lwbnewbuf(  omsg  )  )  return!  0  ); 
if  (  ! lwbgetdatbuf (  omsg,  Aodata  )  )  return!  0  ); 

/*  send  message  on  its  way  and  release  the  input  */ 
lwbwrite!  omsg,  out  ); 
lwbrelease!  in  ); 
return!  1  ) ; 
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} 

/*************************  ****************************************  ^ 

ufrun(  in,  out  ) 
int  in,  out; 

{ 

register  int  i,  j,  k; 

/*  sanity  check  */ 

if  (  !in  !!  !omsg  )  return (  0  ); 

if  (  1 lwbgetparam(  omsg,  0,  LWBMSGLINKS,  Alinks  )  )  retum(  0  ); 
if  (  links  >  0  ) 

{ 

lwbunget(  in  ) ; 
return (  1  ) ; 

} 

if  (  ! lwbgetdatbuf (  in,  Aidata  )  )  retum(  0  ); 

/****  PROCESS  THE  DATA  ***/ 

for  (  i  =  0;  i  <  2 fc  i++)  /*  process  first  20  data  points  in  buffer 

V 

{ 

otemp[i]=0.0;  /*  initialize  otempfi]  */ 

for(j=0;  j  <  i+1;  j-H-) 

otemp[i]  =  otemp[i]+B[ j]  *  (double)  idata[i-j] ; 
for  ( j=i+l ;  j  <21;  j-H-) 
otempfi]  =  otemp[i]+B[ j]  *  stor [20-j+i] ; 
odatafi]  =  (float) (otempfi] ) ; 

} 

for  (  i  =  20;  i  <  count;  i++)  /*  process  remainder  of  buffer  */ 

{ 

otemp[i]=0.0;  /*  initialize  otemp[i]  */ 

for( j=0;  j  <  21;  j++) 

otemp[i]  =  otempfi] +B [j]* (double) idata[i-j] ; 
odatafi]  =  (float) (otempfi] ) ; 

} 

for  (  i  =  0;  i  <  20;  i++)  /*  store  values  for  next  buffer  */ 

storfi]  =  (double) idata[count-20+i] ; 
lwbwrite(  omsg,  out  ); 
lwbrelease(  in  ) ; 
return (  1  ) ; 

} 

j* **********************  STOP  **«*******»*******»*********»*******/ 

ufstop(  in,  out  ) 
int  in,  out; 

{ 

if  (  !in  !!  lout  )  return (  0  ); 
lwbwrite(  omsg,  out  ); 
lwbrelease(  in  ); 
return (  1  ) ; 

> 

j ***********************  ***********************************/ 

ufmenu(  key  ) 
int  key; 
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int  i; 

if  (  key  =  1  ) 

{ 

lwbgetint(  "Enter  Sample  Rate(228  or  456)",  228,  456,  Arate) ; 
if  (  ! lwbupdatemenu(  1,  "%d",  rate)  )  retum(  0  ); 

} 

if  (rate=228) 

for(i=0;  i  <  21;  i++)  B [i ]  =  b_228[i]; 
if  (rate=456) 

for(i=0;  i  <  21;  i++)  B[i]  =  b_456[i]; 
return (  1  ) ; 

} 

/***************#**********  DROP  *********************************/ 

ufdropO 

{ 

if  (  !lwbfreemsg(  Aomsg  )  )  retum(  0  ); 
return (  1  ) ; 

} 
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/************** *********************************************************** 

LOP ASS. c 

This  is  an  implementation  of  a  20  order  lopass  filter  whose  coefficients 
are  produced  by  matlab  "firl"  function. 

***************»*****«»***************«***»******»**********»*********/ 

static  char  SccsID[]  =  "@(#)uf test .c8. 1  (MASSOOMP)  2/5/90"; 

#include  "lwbufunc.h" 

#include  "errors. h" 

/*  input  and  output  data  pointers  */ 
float  ‘idata,  ‘odata; 

int  omsg=0f  links=0,  states,  itype=0,  otype=LWBFLOAT ,  dimensions, 
sizes [3]={0, 0,0} ,  counts,  rate=228,  dxS; 
double  b_228 [21]  =  {-5.51905550574465e-04,  7. 85822093 157633e-04, 

4 . 1 8562857927650e-03 ,  1.1 8776604 1 47200e-02 , 

2 . 54520047550 1 6 1  e-02 ,  4.511 29392503578e-02 , 

6 . 929225487645 1 4e-02 ,  9 . 47857 1 0 1 643358e-02 , 

1 . 174131 8500 1 872e-0 1 ,  1 . 33031630503429e-01 , 

1 . 38608237639420e-0 1 ,  1 . 3303 1 630503429e-0 1 , 

1 . 1741 31 85001 872e-01 ,  9. 47857101 643358e-02, 

6 . 929225487645 1 4e-02 ,  4.511 29392503578e-02 , 

2 . 54520047550 1 6 1 e-02 ,  1.1 8776604 1 47200e-02 , 

4 . 1 8562857927650e-03 ,  7 . 85822093 1 57633e-04 , 

-5.51 905550574465e-04} 

b_456[21]  =  {  5 . 02260 1 74274603e-03 ,  7. 151 15484736450e-03, 

1.281 3237 1 720883e-02 ,  2 . 22206445774868e-02 , 

3.49701 1 95466541 e-02,  5. 00555223943375e-02, 

6 . 59856590758503e-02 ,  8.0991319841 7360e-02 , 

9 . 32883970 1 05506e-02 ,  1.01 352993037039e-0 1 , 

1.04161 357873986e-0 1 ,  1.01 352993037039e-0 1 , 

9.328839701 05506e-02 ,  8.0991319841 7360e-02 , 

6 . 59856590758504e-02 ,  5 . 00555223943375e-02 , 

3 . 4970 1 1 9546654 1 e-02 ,  2 . 22206445774868e-02 , 

1.281 3237 1 720883e-02 ,  7.1511 5484736450e-03 , 

5 . 02260 1 74274603e-03} 

stor[20]={l, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0,0} , 

B[21],  otemp[500]; 

/******»*****************************«***»*********«********************* 

*  »  Entry  point  from  LWB  —  must  be  first  function  defined  in  file  <<  * 
************************************************************************ / 
uftest(  action,  inA,  inB,  inC,  inD,  outX,  outY,  debug  ) 
int  ‘action,  *inA,  *inB,  *inC,  *inD,  outX,  outY,  ‘debug; 

{ 

int  status; 
dx  =  ‘debug; 
switch  (  ‘action  ) 

{ 

case  SETUPCMD: status  =  ufsetup(  *inA,  outX  );  break; 
case  RUNCMD:  status  =  ufrun(  *inA,  outX  );  break; 
case  STOPCMD: status  =  ufstop(  ‘inA,  outX  );  break; 
case  MENUCMD: status  =  ufmenu(  *inA  );  break; 
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case  DROPCMD: status  =  ufdropO ;  break; 
default:  break; 

} 

retum(  status  ); 

} 

uf setup (  in,  out  ) 
int  in,  out; 

{ 

int  i ; 

/*  set  up  menu  label  and  initial  value  */ 

if  (  ! lwblabelmenu(  1,  "Sample  Rate"  )  )  retum(  0  ); 
if  (  ! lwbupdatemenu(  1,  "%d",  rate)  )  retum(  0  ); 
if  (rate=228) 

for(i=0;  i  <  21;  i++)  B[i]  =  b_228[i]; 
if  (rate=456) 

for(i=0;  i  <  21;  i++)  B[i]  =  b_456[i]; 

/*  get  dimensionality  and  size  of  input  message  data  */ 
lwbgetparam(  in,  0,  LWBMSGDIM,  Adimension  ); 
if  (  dx  )  printf(  "dimension=9ad\n",  dimension  ); 
lwbgetparam(  in,  0,  LWBMSGSIZES,  sizes  ); 

if  (  dx  )  printf(  "sizes=%d,%d,%d\n" ,  sizes [0],  sizes [1],  sizes [2]  ); 
if  (  dimension  =  2  )  count  =  sizes [0]  *  sizes [1]; 

else  if  (  dimension  =  1  )count  =  sizes [0]; 
if  (  dx  )  printf(  "count=$d\n",  count  ); 

/*  check  input  data  type  */ 

lwbgetparam(  in,  0,  LWBMSGTYPE,  Ai type  ); 
if  (  dx  )  printf(  " i type=^d\n" ,  itype  ); 
if  (  itype  !=  5  ) 

{ 

lwberror(  ERROR,  "Input  message  data  type  must  be  float"  ); 
return (  0  ) ; 

} 

/*  create  a  new  message  for  output  data  */ 
if  (  !lwbnewmsg(  Aomsg  )  )  retum(  0  ); 

/*  initialize  message  */ 

if  (  !lwbcopymsg(  in,  omsg  )  )  retum(  0  ); 
if  (  ! lwbsamerange(  in,  omsg  )  )  retum(  0  ); 
if  (  ! lwbsamedomain(  in,  omsg  )  )  retum(  0  ); 

/*  create  a  data  buffer  for  output  message  */ 
if  (  !lwbnewbuf(  omsg  )  )  retum(  0  ); 
if  (  !  lwbgetdatbuf  (  omsg,  Aodata  )  )  retum(  0  ); 

/*  send  message  on  its  way  and  release  the  input  */ 
lwbwrite(  omsg,  out  ); 
lwbrelease(  in  ); 
return (  1  ) ; 

> 

/**«*********«************  run  ****************************************/ 

ufrun(  in,  out  ) 
int  in,  out; 

{ 

register  int  i,  j,  k; 
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/*  sanity  check  */ 

if  C  !in  ,'!  Jomsg  )  retum(  0  ); 

if  (  !lwbgetparam(  omsg,  0,  LWBMSGLINKS,  Alinks  )  )  retum(  0  ); 
if  (  links  >  0  ) 

{ 

lwbunget(  in  ) ; 
return (  1  ); 

} 

if  (  ! lwbgetdatbuf (  in,  &idata  )  )  retum(  0  ); 

/****  PROCESS  THE  DATA  ***/ 

for  (  i  =  0;  i  <  20;  i++)  /*  process  first  20  data  points  in  buffer 

V 

{ 

otemp[ij=0.0;  /*  initialize  otemp[ij  */ 

for( j=0;  j  <  i+1;  j++) 

otempfi]  =  otemp[i]+B{j]  *  (double) idata[i-j] ; 
for(j=i+l;  j  <  21;  j++) 

o temp[i]  =  otemp[i]+B[j]  4  stor[2 0-j+i] ; 
odata[i]  =  (float) (otemp(i]> ; 

> 

for  (  i  =  20;  i  <  count;  i++)  /*  process  remainder  of  buffer  */ 

{ 

otemp[i]=0.0;  /*  initialize  otemp[i]  */ 

for( j=0;  j  <  21;  j++) 

otemp[i]  =  otemp[i]+B[j] 4 (double) idata[i-j] ; 
odata[i)  =  (float) (otemp[i)) ; 

} 

for  (  i  =  0;  i  <  20;  i++)  /*  store  values  for  next  buffer  */ 

stor[i]  =  (double) idata[count-20+i] ; 
lwbwrite(  omsg,  out  ); 
lwbrelease(  in  ); 
return (  1  ); 

} 

/*************#*********  STOP  ************************************  j 

ufstop(  in,  out  ) 
int  in,  out; 

( 

if  (  Jin  II  lout  )  return (  0  ); 
lwbwrite(  omsg,  out  ); 
lwbrelease(  in  ); 
re turn (  1  ) ; 

} 

/**«**********«*********  menu  ***********************************/ 

ufmenu(  key  ) 
int  key; 

{ 

int  i ; 

if  (  key  ==  1  ) 

{ 

lwbgetint(  "Enter  Sample  Rate(228  or  456)",  228,  456,  &rate) ; 
if  (  ! lwbupdatemenu(  1,  "%d",  rate)  )  return(  0  ); 
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} 

if  (rate=228) 

for(i=0;  i  <  21;  i++)  B[i]  =  b_228[i]; 
if  (rate=456) 

for(i=0;  i  <  21;  i++)  B[i]  =  b_456[i]; 
return (  1  ) ; 

} 

!% *************************  DROP  *********************************/ 
ufdropO 
{ 

if  (  !  lwbfreemsg(  &omsg  )  )  retum(  0  ); 
retum(  1  ) ; 

} 
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/ft********************************************************************* 

FDEMODl.c 

This  program  forms  a  user  function  in  Lab  Workbench  to  demodulate  the 
input  data  stream  and  send  two  output  data  streams  at  x  and  y  consisting 
to  the  real  and  imaginary  parts,  respectively.  The  user  can  select  the 
carrier  and  sample  frequency  via  menu.  Note:  The  data  MUST  be 

demodulated  at  exactly  the  frequency  of  the  received  carrier  frequency  for 
the  sequence  removal  to  function  optimally. 

static  char  SccsID[]  =  "8(#)uf test.c8. 1  (MASSCOMP)  2/5/90"; 

#include  "lwbufunc.h" 

#include  "errors. h" 

#include  <math.h> 

/*  input  and  output  data  pointers  */ 
float  ‘idata,  ‘odataX,  ‘odataY; 

int  omsgX=0,  omsgY=0,  linksX=0,  linksY=0,  itype=0,  idimension=0, 
odimension=l , 

insizes [3]={0, 0,0} ,  outsizes[3],  count=0,dx=0; 
static  n=0; 

double  fc=56.983,  fs=228.0,  theta=0.0,  pi,  doppler=0.0,  doppler_bin=0.0, 
sndve 1=1 500.0; 

/*  ********** ************************************************************* 

*  »  Entry  point  from  LWB  —  must  be  first  function  defined  in  file  <<  * 

************************************************************************  / 

uftest(  action,  inA,  inB,  inC,  inD,  outX,  outY,  debug  ) 
int  ‘action,  *inA,  *inB,  *inC,  *inD,  outX,  outY,  ‘debug; 

{ 

int  status; 
dx  =  ‘debug; 
switch  (  ‘action  ) 

{ 

case  SETUPCMD: status  =  ufsetup(  *inA,  outX,  outY  );  break; 
case  RUNCMD:  status  =  ufrun(  ‘inA,  outX,  outY  );  break; 
case  STOPCMD: status  =  ufstop(  *inA,  outX,  outY  );  break; 
case  MENUCMD: status  =  ufmenu(  *inA  );  break; 
case  DROPCMD: status  =  ufdropO ;  break; 
default:  break; 

> 

re turn (  status  ) ; 

} 

/% *****************************  SETUP  *******************************/ 

uf setup (  in,  outX,  outY  ) 
int  in,  outX,  outY; 

{ 

double  asin(); 
n=0; 

odimension=l ; 

pi  =  2  *  asin(l .0) ; 

/*  set  up  menu  label  and  initial  value  */ 

if  (  ! lwblabelmenu(  1,  "Carrier  Frequency"  )  )  retum(  0  ); 
if  (  ! lwbupdatemenu(  1,  "%f",  fc  )  )  retum(  0  ); 
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from  setup  =%f\n". 


if  (  !  lwb  label  menu  (  2,  "Sample  Frequency"  )  )  retum(  0  ); 
if  (  ! lwbupdatemenu(  2,  "%f",  fs  )  )  retum(  0  ); 
if  (  !  1 wb label menu (  3,  "Doppler  (m/s)"  )  )  retum(  0  ); 
if  (  ! lwbupdatemenu(  3,  "%f",  doppler  )  )  retum(  0  ); 
if  (  ! lwblabelmenu(  4,  "Sound  Speed  (m/s)"  )  )  retum(  0  ) 
if  (  !  lwbupdatemenu(  4,  "%f",  snd_vel  )  )  retum(  0  ); 

/*  compute  theta  */ 

theta  =  2  *  pi  *  (  fc  /  fs  ); 

if  (  dx  )  printf(  "demod  theta  from  setup  =%f\n",  theta  ); 
/*  compute  doppler  correction  */ 

doppler_bin  =  pi  *  doppler  /  (2  *  snd_vel); 
if  (  dx  )  printf(  "demod  doppler_bin 
doppler_bin  ) ; 

/*  check  input  messages  */ 

/*  get  dimensionality  of  input  message  data  */ 
lwbgetparam(  in,  0,  LWBMSGDIM,  ftidimension  ); 
if  (  dx  )  printf(  "demod  idimension3^^" ,  idimension  ); 

/*  get  size  of  input  message  data  */ 
lwbgetparam(  in,  0,  LWBMSGSIZES,  insizes  ); 
count  =  insizes [0] ; 

if  (  dx  )  printf(  "demod  count=55d\n" ,  count  ); 

/*  check  input  data  type  */ 
lwbgetparam(  in,  0,  LWBMSGTYPE,  ftitype  ); 
if  (  dx  )  printf(  "demod  itype=56d\n" ,  itype  ); 

(  itype  !=  5  ) 


if 

{ 


lwberror(  ERROR,  "Input  message  data  type  must  be  FLOAT"  ); 
re turn (  0  ) ; 


} 

/*  create  a  new  message  for  output  data  */ 

if  (  !lwbnewmsg(  ftomsgX  )  )  return(  0  ); 

if  (  ’lwbnewmsg(  ftomsgY  )  )  retum(  0  ); 

/*  copy  msg  header  info  from  input  to  output  message  */ 
if  (  !lwbcopymsg(  in,  omsgX  )  )  retum(  0  ); 

if  (  !lwbcopymsg(  in,  omsgY  )  )  return(  0  ); 
outsizes [0]=insizes [0] ; 
outsizes [lj=l ; 
outsizes [2] =0; 

if  (  !lwbsetparam(  omsgX,  0,  LWBMSGSIZES,  outsizes  )  )  retum(  0  ); 

if  (  ! lwbsetparam(  omsgY,  0,  LWBMSGSIZES,  outsizes  )  )  retum(  0  ); 

if  (  ! lwbsetparam(  omsgX,  0,  LWBMSGDIM,  Aodimension  )  )  retum(  0  ); 

if  (  ! lwbsetparam(  omsgY,  0,  LWBMSGDIM,  Aodimension  )  )  return(  0  ); 

/*  set  output  range  parameters  */ 
if  (  ! lwbnewrange(  omsgX,  1  )  )  retum(  0  ); 

if  (  ! lwbnewrange (  omsgY,  1  )  )  retum(  0  ); 

if  (  !  lwbcopyrange(  in,  omsgX  )  )  retum(  0  ); 
if  (  !lwbcopyrange(  in,  omsgY  )  )  retum(  0  ); 
if  (  !lwbnewdomain(  omsgX  )  )  return!  0  ); 

if  (  ! lwbnewdomain(  omsgY  )  )  retum(  0  ); 

if  (  ! lwbcopydomain(  in,  omsgX  )  )  retum(  0  ); 
if  (  ! lwbcopydoraain(  in,  omsgY  )  )  return!  0  ); 
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/*  create  a  data  buffer  for  output  message  */ 
if  (  !lwbnewbuf(  omsgX  )  )  retum(  0  ); 

if  (  !lwbnewbuf(  omsgY  )  )  retum(  0  ); 
if  (  ! lwbgetdatbuf (  omsgX,  AodataX  )  )  return(  0  ); 

if  (  ! lwbgetdatbuf (  omsgY,  AodataY  )  )  retum(  0  ); 
lwbgetparam(  omsgX,  0,  LWBMSGDIM,  &odi mens  ion  ); 
if  (  dx  )  printf(  "demod  odimension=%d\n" ,  odimension  ); 

/*  send  message  on  its  way  */ 
lwbwrite(  omsgX,  outX  ); 
lwbwrite(  omsgY,  outY  ); 

/*  and  release  the  input  */ 
lwbrelease(  in  ); 
return(  1  ) ; 

} 

/**************************  **********************************  / 

ufrun(  in,  outX,  outY  ) 
int  in,  outX,  outY; 

{ 

register  int  i,  j; 
static  m; 

double  temp,  cos(),  sin(); 

/*  sanity  check  */ 

if  (  !in  ij  !omsgX  ||  lomsgY  )  return(  0  ); 

if  (  ! lwbgetparam(  omsgX,  0,  LWBMSGLINKS,  AlinksX  )  )  return (  0  ); 

if  (  ! lwbgetparam(  omsgY,  0,  LWBMSGLINKS,  filinksY  )  )  retum(  0  ); 
if  (  (linksX  >  0)  | |  (linksY  >  0)  ) 

{ 

lwbunget(  in  ) ; 
return (  1  ) ; 

} 

if  (  ! lwbgetdatbuf (  in,  Aidata  )  )  retum(  0  ); 
for  (i=0;  i  <  count;  i++) 

{ 

temp  =  (double) idata[i] ; 
odataX[i]  =  (float) (temp  *  cos (n* theta) ) ; 
odataY[i]  =  (float) (temp  *  sin(n*theta)) ; 
n=n+l ; 

} 

m++; 

lwbwrite(  omsgX,  outX  ); 

lwbwrite(  omsgY,  outY  ); 
lwbrelease(  in  ); 
return (  1  ) ; 

} 

/*****************************  STOP  *************************«*******/ 

ufstop(  in,  outX,  outY  ) 
int  in,  outX,  outY; 

{ 

if  (  !in  i!  JoutX  ',!  !outY  )  retum(  0  ); 
lwbwrite(  omsgX,  outX  ); 
lwbwrite(  omsgY,  outY  ); 
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} 


lwbrelease(  in  ); 
return (  1  ) ; 


y* ***************************  **********************  ***********/ 

ufmenu(  key  ) 
int  key; 

{ 

if  (  key  —  1  ) 

{ 

lwbgetdouble(  "Enter  Carrier  Frequency",  1.0,  1000.0,  &fc  ); 
if  (  ! lwbupdatemenu(  1,  "%f",  (float)fc  )  )  retum(  0  ); 

} 

if  (  key  =  2  ) 

{ 

lwbgetdouble(  "Enter  Sample  Frequency",  1.0,  1000.0,  &fs 

) ; 

if  (  ! lwbupdatemenu(  2,  "%f",  (float)fs  )  )  retum(  0  ); 

} 

if  (  key  =  3  ) 

{ 

lwbgetdouble(  "Enter  Total  Doppler  (m/s) ",-5.0, 5.0, 

&doppler) ; 

if  (  ! lwbupdatemenu(  3,  "%f",  (f loat)doppler)  )  retum(  0 

); 

} 

if  (  key  =  4  ) 

{ 

1 wbg e t doub 1 e ( " En t e r  Sound  Velocity 
(m/s) " , 1400.0, 1500.0, &snd_vel) ; 

if  (  ! lwbupdatemenu(  4,  "%f",  (f loat)snd_vel)  )  retum(  0  ); 

} 

theta  =2*pi*(fc/fs); 
doppler_bin  =  pi  *  doppler  /  (2  *  snd_vel) ; 
return (  1  ) ; 

} 

/***« ***************************  drop  *******************************/ 

ufdropQ 

{ 

if  (  !lwbfreemsg(  AomsgX  )  )  return(  0  ); 

if  (  !lwbfreemsg(  AomsgY  )  )  retum(  0  ); 
return (  1  ) ; 

} 
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ft  ************************  *******  *******  ******************************** 


FMETl.c 

This  is  the  sequence  removal  program  which  calls  a  subroutine  developed  by 
Metzger  to  do  the  pulse  compression.  This  program  forms  a  user  function 
in  Lab  Workbench  and  allows  menu  selection  of  the  sample  freq.  and  law  of 
the  code.  Note  that  LWB  doesn't  let  you  enter  octal  so  enter  the  decimal 
law  of  the  code  and  it  will  be  converted  and  displayed  as  the  octal  law  on 
the  menu. 

*****«*«4 


static  char  SccsID[]  =  "@(#)uftest .c8. 1  (MASSCOMP)  2/5/90”; 
#include  "lwbufunc.h" 

#include  "errors. h" 

#include  <malloc.h> 

#include  <stdio.h> 

#include  <math.h> 

#def ine  PI  3.1415926536 


/*  input  and  output  data  pointers  */ 
float  *indataA,  ♦indataB,  ‘outdata; 

double  iangle=45.0,  angle,  factor,  fact_i; 


float 


slope, 
baseline, 
df  s=456 , 
df  c=57 , 
cycles_dig=5, 
thresh_fact=50; 


/*  sample  freq  (changable  by  menu)  */ 

/*  carrier  freq  */ 

/*  carrier  cycles  per  digit  */ 

/*  printout  if  output=thresh_fact*out_ave  */ 


/*  (global)  variables 
int  law=351, 
initial=l , 
counts, 
dem_dig=20, 
seq_len=l , 
degree , 
outmsg=0 , 
links=0, 
intypeA=0, 
intypeB=0, 
indimA=0, 
indimB=0, 
insizesA(3] , 
insizesB[3] , 
osizes[3] , 
gain, 

conv_gain, 
step=l , 
dx=0, 

scram [2047] , 
unscram [2047] , 
dec_data [83000] , 


(put  these  in  a  macrofile  later)  */ 

/*  decimal  law.. program  converts  to  octal  */ 
/*  initial  register  load  */ 

/*  data  points  in  one  input  buffer  */ 

/*  data  points  in  one  phase  digit  */ 

/*  data  points  in  one  m-sequence .. computed  */ 
/*  number  of  registers  in  sequence  generator  */ 


/*  LWB  inputA  data  type  (ie  float=5)  */ 

/*  dimension  of  inputA  buffer  array  */ 

/*  array  to  hold  input  buffer  size  data  */ 


/*  value  of  least  sig.  bit  of  data  */ 

/*  decimation  in  time  factor  if  used  */ 

/*  flag  to  execute  debug  statements  */ 

/*  array  to  hold  scrambled  input  data  */ 

/*  array  to  hold  unscrambled  input  data  */ 
/*  array  to  accumulate  input  data  */ 
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seq_sum[80] , 


cbuf [4100] ; 

static  long  k=0, 

1=0, 
mark_l , 
num_wrap; 


/*  array  to  hold  the  0th  item  of  each  sequence 
passed  to  hadamard  routine  so  bias  can  be 
accumulated  */ 

/*  array  to  temporarily  hold  one  seq_len  of 

data  to  be  passed  to  hadamard  routine  */ 
/*  data  point  counter  */ 

/*  buffer  counter  */ 

/*  marks  first  excess  data  in  k  array  */ 

/*  number  of  exess  data  points  in  k  array  */ 


/*********************************************************««*******««*««* 

*  >>  Entry  point  from  LWB  —  must  be  first  function  defined  in  file  <<  * 

************************************************************************ i 


uftest(  action,  inA,  inB,  inC,  inD,  outX,  outY,  debug  ) 
int  ‘action,  *inA,  *inB,  *inC,  *inD,  outX,  outY,  ‘debug; 

{ 

int  status; 
dx  =  ‘debug; 
switch  (  ‘action  ) 

{ 

case  SETUPCMD: status  =  ufsetup(  ‘inA,  *inB,  outX  );  break; 
case  RUNCMD:  status  =  ufrun(  *inA,  *inB,  outX  );  break; 
case  STOPCMD: status  =  ufstop(  *inA,  *inB,  outX  );  break; 
case  MENUCMD: status  =  ufmenu(  *inA,  *inB  );  break; 
case  DROPCMD:  status  =  ufdropQ ;  breed;; 
default:  break; 

} 

return (  status  ) ; 


!% *******************************  SETUP  ft**************************/ 

ufsetup(  inA,  inB,  out  ) 
int  inA,  inB,  out; 


{ 

int 

unsigned 
float 
double 
initial=l ; 


temp=law,  rev_law=0,  end_bit,  i,  j; 
rev_initial,  ss_contents,  ms_contents; 
chanhilim; 

denomin,  taneing,  tansqr; 


mark_l=0;  /*  initialize  excessdata  marker  •/ 

num_wrap=0;  /*  initialize  excess  data  counter  ♦/ 
k=0;  /*  initialize  data  counter  */ 

1=0;  /*  initialize  buffer  counter  ♦/ 

/*  menu  label  and  initial  value  (note  need  351  to  get  correct  537)  */ 
if  (  !lwblabelmenu(  1,  "OCTAL  LAW"  )  )  retum(  0  ); 
if  (  ! lwbupdatemenu(  1,  "%o",  law  )  )  retum(  0  ); 
if  (  !lwblabelmenu(  2,  "Initial  Load"  )  )  retum(  0  ); 
if  (  ! lwbupdatemenu(  2,  "%d",  initial  )  )  retum(  0  ); 
if  (  !  1  wb  label  menu  (  3,  "Phase  Angle"  )  )  retum(  0  ); 
if  (  ! lwbupdatemenu(  3,  "%f",  iangle  )  )  return(  0  ); 
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if  (  ! lwblabelmenu(  4,  "Sample  Freq"  )  )  retum(  0  ); 
if  (  ! lwbupdatemenu(  4,  "%f",  dfs)  )  retum(  0  ); 

/*  compute  demodulates  per  digit  */ 
dem_dig  =  (int)((dfs  /  dfc) *cycles_dig) ; 
if  (  !  1  wb  label  menu  (  5,  "Threshold"  )  )  retum(  0  ); 
if  (  ! lwbupdatemenu(  5,  "%f",  thresh_fact)  )  return(  0  ); 

/*  check  input  data  */ 

/*  get  dimensionality  and  size  of  input  message  data  */ 

lwbgetparam(  inA,  0,  LWBMSGDIM,  AindimA  ); 

lwbgetparam!  inB,  0,  LWBMSGDIM,  AindimB  ); 

lwbgetparam!  inA,  0,  LWBMSGSIZES,  insizesA  ); 

lwbgetparam!  inB,  0,  LWBMSGSIZES,  insizesB  ); 

if  (  indimA  =  2  )  counts  =  insizesA[0]  *  insizesAfl]; 

else  if  (  indimA  =  1  )counts  =  insizesA[0]; 

/*  check  input  data  type  */ 

lwbgetparam!  inA,  0,  LWBMSGTYPE,  AintypeA  ); 

if  (  intypeA  !=  5  ) 

{ 

lwberror(  ERROR,  "Input  A  message  data  type  must  be  FLOAT"  ); 
return!  0  ) ; 

} 

lwbgetparam!  inB,  0,  LWBMSGTYPE,  AintypeB  ); 
if  (  intypeB  !=  5  ) 

{ 

lwberror!  ERROR,  "Input  B  message  data  type  must  be  FLOAT"  ); 
return!  0  ) ; 

} 

/*  get  gain,  convgain,  slope,  and  baseline  data  from  msg  header  */ 

lwbgetparam!  inA,  0,  LWBCHANGAIN,  Again  ); 

lwbgetparam!  inA,  0,  LWBCHANCONVGA I N ,  Aconvgain  ); 

lwbgetparam!  inA,  0,  LWBCHANSLOPE ,  Aslope  ); 

lwbgetparam!  inA,  0,  LWBCHANBASELINE,  Abaseline  ); 

/*  compute  sequence  length  (digitz/seq.)  and  polynomial  degree  */ 
temp=law; 
seqL_len=l ; 
rev_initial  =  0; 
degree  =  0; 

while  (temp»=l) 

{ 

seqL_len  <<=  1 ; 

degree-H- ; 

> 

/*  compute  scram  and  unscram  arrays  */ 
end_bit=seq_len — »1 ; 
for  (i=0;  i  <  degree;  i++) 

{ 

revinitial  =  (rev_initial  <<  1)  !  (initialAl); 
initial  >>=  1; 

} 

ms_contents  =  1 ; 
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ss_contents  =  rev_initial; 
for  (i  =  0;  i  <  seq_len;  i++) 

{ 

scram [i]  =  ss_contents; 
unscram[i]  =  ms_contents; 
temp  =  ss_contents  &  law; 
ss_contents  >>=  1 ; 

do  {if  (temp  &1)  ss_contents'=end_bit;}  while  (temp>>=l); 
if  (ms_contents  &1) 

ms_contents  =  (ms_contents~law)  »  1; 


% 


else 


ms_contents  >>=  1; 

} 

/*  compute  dc  bias  removal  factors  to  be  used  in  run  */ 
angle  =  iangle  *  PI  /  180.0; 
tanang  =  tan(angle) ; 
tansqr  =  tanang  *  tanang; 

denomin  =  ( (double) seq_len  *  (double) seq_len  +  tansqr); 
fact_r  =  (tansqr  -  (double) seq_len)  /  denomin; 
fact_i  =  ( (double) seq_len  +1.0)  *  tanang  /  denomin; 

/*  prepare  output  message  */ 

/*  copy  input  general  info  to  output  message  */ 
if  (  !lwbnewmsg(  fioutmsg  )  )  retum(  0  ); 
if  (  !lwbcopymsg(  inA,  outmsg  )  )  retum(  0  ); 

/*  change  output  data  size  */ 
osizes [0]=  seq_len  *  dem_dig; 
osizes [1]=  1; 
osizes [2]=  0; 

lwbsetparam(  outmsg,  0,  LWBMSGSIZES,  osizes  ); 

/*  copy  input  message  range  info  to  output  message  */ 
if  (  ! lwbnewrange(  outmsg,  1  )  )  retum(  0  ); 
if  (  ! lwbcopy range (  inA,  outmsg  )  )  return(  0  ); 

/*  set  new  domain  parameters  */ 

if  (  ! lwbnewdomain(  outmsg  )  )  retum(  0  ); 

if  (  ! lwbcopydomain(  inA,  outmsg  )  )  retum(  0  ); 

if  ( ! lwbgetparam(  outmsg,  -1,  LWBCHANHILIM,  &chanhilim))  retum(0); 

chanhilim  *=  seq__len; 

if  (!lwbsetparam(  outmsg,  -1,  LWBCHANHILIM,  Achanhilim))  return (0); 
lwbgetparam(  outmsg,  0,  LWBMSGSIZES,  osizes  ); 

/*  create  a  data  buffer  for  output  message  */ 

if  (  !lwbnewbuf(  outmsg  )  )  return!  0  ); 

if  (  ! lwbgetdatbuf (  outmsg,  &outdata  )  )  return!  0  ); 

/*  send  message  on  its  way  and  release  the  input  */ 

lwbwrite(  outmsg,  out  ); 

lwbrelease(  inA  ); 

lwbrelease(  inB  ); 

return(  1  ) ; 

} 

/************»***********************  RUN  ****************************/ 

ufrun(  inA,  inB,  out  ) 
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int  inA,  inB,  out; 

{ 

int  m,  /*  index  counter  for  stp  */ 

hadamard() ,  /*  hadamard  transform  routine  */ 

*scram_ptr,  /*  pointer  to  scrambled  data  array  */ 

*unscram_ptr,  /*  pointer  to  unscrambled  data  array  */ 

stp,  /*  data  items/digit  since  must  do  hadamard 

correlation  once  for  every  data  item 
(ie  20*2  times  for  re  &  im  parts  @  228 

hz)  */ 

i,  /*  general  counter  index  */ 

j,  /*  general  counter  index  */ 

cor_r,  /*  real  correction  factor  */ 

cor_i,  /*  imag  correction  factor  */ 

*ptr,  /*  ptr  to  main  data  array  (dec_data)  */ 

*data_ptr,  /*  ptr  to  current  data  item  of  digit 

(ie  lst.or  2nd,  etc)  */ 

*move_ptr,  /*  ptr  to  temp  array  of  seq_len  digits  to  be 

passed  to  hadamard  routine  (ie  all  the 

1st, 

or  2nd,  etc  data  items  of  each  digit  */ 
*sum_ptr;  /*  ptr  to  seq_sum  array  which  holds  the  0th 

item  of  each  sequence  passed  to  hadamard 
to  accumulate  bias  */ 

double  xt,  yt,  sqrt(),  r_data,  i_data; 

float  out_sum,  outave,  pts_sum; 

/*  sanity  check  */ 

if  (  !inA  !!  loutmsg  )  retum(  0  ); 

if  (  !inB  ||  loutmsg  )  retum(  0  ); 

if  (  ! lwbgetpar£im(  outmsg,  0,  LWBMSGLINKS,  ftlinks  )  )  return(  0  ); 
if  (  links  >  0  ) 

{ 

lwbunget(  inA  ) ; 
lwbunget(  inB  ); 
return (  1  ) ; 

} 

if  (  ! lwbgetdatbuf (  inA,  AindataA  )  )  retum(  0  ); 

if  (  ! lwbgetdatbuf (  inB,  ftindataB  )  )  return(  0  ); 

/*  enter  demodulates  into  dec_data  array  one  buffer  at  a  time 

(r0,i0,rl , il . . r 1 9 , i 1 9)  until  a  full  sequence  is  recieved 

(dem_dig*seq_len*2  =  10200  data  items  for  law537  @  228hz) 

(r0,i0, . . .rl0199, i 10199)  */ 

if  (k  <  (seq_len  *  dem_dig  *  2  /step))  /*  ie  if  less  than  a  full 
sequence 

has  been  recieved  */ 

( 

if  ((num_wrap  >0)  &&  (k=0))  /*  if  this  is  1st  buffer  of  a  new  seq 

and  there  was  a  remainder  from  the 
last  sequence  */ 

{ 

for  (  i  =  0;  i  <  num_wrap;  i  +=  step) 
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{ 


dec_data [k++j  =  dec_data[mark_l++] ; 
dec_data  [k-H-]  =  dec_data[mark_l++] ; 

> 

} 

for  (  i  =  0;  i  <  counts;  i  +=  step) 

{ 

dec_data [k++]  =  (long) (le7*indataA[i]) ;  /*  note:  mult  by  number 

*  to  get  large  enough  to 

convert  to  long  int  */ 

dec_data[k++]  =  (long) (le7*indataB[i]) ; 

*  if (k  =  (seq_len  *  dem_dig  *  2  /  step)) 

{ 

mark_l  =  k;  /*  set  marker  for  first  data  of  next  sequence  */ 
num_wrap=((counts-i)-l) ;  /*  number  of  data  to  wrap  to  next  seq*/ 

} 

} 

} 

if  (k  >=  (seq_len  *  dem_dig  *  2  /step))  /*  after  a  full  sequence 

is  received...  */ 

{ 

ptr  =  decdata;  /*  ptr  to  main  data  array  (dec_data)  */ 

sum_ptr  =  seq_sum;  /*  assign  sum_ptr  to  seq_sum  array  */ 

move_ptr  =  cbuf;  /*  assign  move_ptr  to  cbuf  array  */ 

stp  =  dem_dig  *  2;  /*  define  stp  as  data  items  per  digit 

(real  &  im)  */ 

/*  Take  the  nth  data  element  of  each  digit  in  the  m-sequence , 
scramble  them,  perform  hadamard  transform,  and  then  unscramble 
them.  Repeat  for  each  successive  element  in  the  digits.  */ 
for  (m  =  0;  m  <  stp;  m++) 

{ 

*move_ptr  =  0;  /*  initial  entry  is  zero  */ 

data_ptr  =  ptr  +  m; 

scram_ptr  =  scram; 

for  (i=0;  i  <  seq_len;  i++) 

{ 

•(movejptr  +  *scramjptr++)  =  *data_ptr;  /*  get  data  values  */ 
data_ptr  +=  stp; 

r  > 

hadamard (degree,  move_ptr) ; 

*sum_ptr++  =  *move_ptr; 

^  data_ptr  =  ptr  +  m; 

unscram_ptr  =  unscram; 
for  (i=0;  i  <  seq_len;  i++) 

{ 

*data_ptr  =  *(move_ptr  +  *unscram_ptr++) ;  /*  put  results  back*/ 
data_ptr  +=  stp; 

} 

} 

for  (i=0;  i  <  dem_dig;  i++)  /*  bias  removal  loop  */ 

{ 
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xt  =  (double)  seq_sum[i+i]  ; 

yt  =  (double) seq_sum [i+i+1] ; 

cor_r  =  (long)(xt  *  fact_r  -  yt  *  fact_i) ; 

cor_i  =  (long)(xt  *  fact_i  +  yt  *  fact_r) ; 

data_ptr  =  ptr  +  i  +  i; 

for  (j  =  0;  j  <  seq_len;  j++) 

{ 

*data_ptr++  -=  cor_r ; 

*data_ptr - -  cor_i; 

data_ptr  +=  (dem_dig  <<  1); 

> 

> 

/*  write  to  output  buffer  */ 
j=0*. 

for  (i=0;  i  <  seq_len  *  dem_dig;  i++) 

{ 

r_data  =  (double) (dec_data[j++]*le~7) ; 
i_data  =  (double) (dec_data[ j++] *le-7) ; 
outdatafi]  =  (float) ((r_data*r_data)+(i_data*i_data)) ; 

/*  if  (  dx  )  printf(  "decoder  outdata[6]=%f\n" ,  outdata[6]) ;*/ 
out_sum  +=  outdatafi]; 

} 

pts_sum  =  (float ) (seq_len  *  dem_dig) ; 
out_ave  =  outsum  /  pts_sum; 

if  a  point  exceeds  thresh_fact  *  out_ave,  print  buffer#,  out_ave 
seconds,  and  value 

for  (i=0;  i  <  seq_len  *  dem_dig;  i++) 

{ 

if  (outdata[i]  >  out_ave  *  thresh_fact) 

{ 

printf ("buffer  #%d  ",1); 

printf ("ave^%4.4e  v.  ",  out_ave) ; 

printf ("%4.3f  sec.  =  %4.4e  v.\n",  i  /  dfs,  outdatafi]); 

} 

} 

/*  write  output  data  message  */ 
lwbwrite(  outmsg,  out  ); 
k  =  0; 

1  =  1+1; 

} 

/*  release  input  message  */ 
lwbrelease(  inA  ); 
lwbrelease(  inB  ); 
return (  1  ) ; 

} 

/  *  HADAMARD  FHT  routine 

*************************************************/ 

#include  <math.h> 
hadamard(n,  data) 
int  n,  *data; 

{ 


4. 
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int  *ptrp,  *ptrl,  *ptr2,  tempi,  temp2; 

unsigned  int  two_step,  parts,  part; 
register  unsigned  int  step,  k; 
step  =  (parts  =  1)  <<  n; 

while  ((step  =  (two_step  =  step)  >>  1)  !=  0) 

{ 

ptrp  =  data; 
part  =  parts; 
do 
{ 

ptrl  =  ptrp; 
ptrp  +=  two_step; 
k  =  step; 
do 
{ 

ptr2  =  ptrl  +  step; 

*ptrl  =  (tempi  =  *ptrl)  +  (temp2  =  *ptr2) ; 

*ptr2  =  tempi  -  temp2; 

ptrl++; 

}  while  ( — k) ; 

}  while  ( — part); 
parts  <<=  1 ; 

} 

return; 

} 


/* ********** ******** *******************  STOP  *********************/ 

ufstop(  inA,  inB,  out  ) 
int  inA,  inB,  out; 

{ 

if  (  linA  !!  lout  )  return (  0  ); 
if  (  linB  {,'  lout  )  retum(  0  ); 
lwbwrite(  outmsg,  out  ); 
lwbrelease(  inA  ); 
lwbrelease(  inB  ); 
return (  1  ) ; 

} 

/* *********************** *************  menu  ***********************/ 

ufmenu(  key  ) 
int  key; 

{ 

if  (  key  =  1  ) 

{ 

lwbget  int  (  "Enter  DECIMAL  Law  ( 1 1 , 351 , 827 , 1 051 , 3047) " ,  1 1 , 3047 , 
&law  ) ; 

if  (  1  lwbupdatemenu(  1,  "%o",  law  )  )  retum(  0  ); 

} 

if  (  key  ==  2  ) 

{ 
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lwbgetint(  "Enter  Initial  Register  Load  (not  zero)",  1,  10000, 
fiinitial  ); 

if  (  ! lwbupdatemenu(  2,  "%d",  initial  )  )  retum(  0  ); 

} 

if  (  key  =  3  ) 

{ 

lwbgetdouble(  "Enter  Phase  Angle",  0.0,  180.0,  Aiangle); 
if  (  ! lwbupdatemenu(  3,  "%f",  iangle  )  )  retum(  0  ); 

} 

if  (  key  =  4  ) 

{ 

lwbgetfloat ("Enter  Sample  Freq",  1.0,  1000.0,  &dfs); 
if  (  !lwbupdatemenu(  4,  "%f",  dfs)  )  retum(  0  ); 

} 

if  (  key  =  5  ) 

{ 

lwbgetfloat ("Enter  threshold  factor",  0.0,  100000.0, 

&thresh_fact) ; 

if  (  ! lwbupdatemenu(  5,  "%f",  thresh_fact  )  )  retum(  0  ); 

} 

return (  1  ); 

} 

/***************»«*******************  DROP  ***********************/ 

ufdropO 

{ 

if  (  !  lwbfreemsg(  Aoutmsg  )  )  retum(  0  ); 
return (  1  ) ; 

} 


/********************************************************************* 

SPEWRIT.C 


This  program  takes  output  from  the  power  spectrum  module  and  finds  the 
highest  peak,  in  a  5  knot  doppler  window.  It  writes  the  peak  frequency  and 
amplitude  (v/Hz)  to  a  file.  It  also  calculates  the  mean  noise  in  a  24  Hz 
band  centered  at  57  Hz  and  writes  it  to  a  file. 

**********************************************************************  f 

static  char  SccsID[]  =  "@(#)uftest.c8. 1  (MASSCOMP)  2/5/90"; 

#include  "lwbufunc.h" 

#include  "errors. h" 

#include  <stdio.h> 

FILE  ‘outfilel,  *outfile2,  *outfile3; 

char  outfilenaml [15]="m"t  outfilenam2[15]="f",  outfilenam3[15]="n" ; 

float  ‘idata,  thresh; 

int  itype=0,  idimension=0,  insizes[3]={0,0,0} ,  count=0,  dx=0; 
static  int  m; 

/************************************************************************ 

*  »  Entry  point  from  LWB  —  must  be  first  function  defined  in  file  <<  * 

************************************************************************ / 

uftest(  action,  inA,  inB,  inC,  inD,  outX,  outY,  debug  ) 
int  ‘action,  *inA,  *inB,  ‘inC,  *inD,  outX,  outY,  ‘debug; 

{ 

int  status; 
dx  =  ‘debug; 
switch  (  ‘action  ) 

{ 

case  SETUPCMD: status  =  ufsetup(  ‘inA  );  break; 
case  RUNCMD:  status  =  ufrun(  *inA  );  break; 
case  STOPCMD: status  =  ufstop(  *inA  );  break; 
case  MENUCMD: status  =  ufmenu!  *inA  );  break; 
case  DROPCMD: status  =  ufdropO ;  break; 
default:  break; 


} 

return (  status  ) ; 


} 

!% *****************************  SETUP  ********************************/ 

ufsetup(  in  ) 
int  in; 

{ 


m=l ; 

/*  set  up  menu  label  and  initial  value  */ 

AAG  filename"  )  )  retum(  0  ); 
"%s",  outfilenaml)  )  retum(  0  ); 
FREQ  filename"  )  )  retum(  0  ); 
"Is",  outfilenam2)  )  return(  0  ); 
'JOISE  filename"  )  )  retum(  0  ); 
"%s",  outfilenam3)  )  retum(  0  ); 
! lwblabelmenu(  4,  "Threshold"  )  )  return!  0  ); 

! lwbupdatemenu(  4,  "%f",  thresh)  )  return!  0  ); 
check  input  messages  */ 

/*  get  dimensionality  of  input  message  data  */ 


if 

( 

! lwblabelmenu!  1, 

if 

( 

! lwbupdatemenu!  1 , 

if 

( 

! lwblabelmenu!  2, 

if 

( 

! lwbupdatemenu!  2, 

if 

( 

! lwblabelmenu!  3, 

if 

( 

! lwbupdatemenu!  3, 

if 

if 
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lwbgetparam(  in,  0,  LWBMSGDIM,  Aidimension  ); 

if  (  dx  )  printf(  "spewrit  idimension=%d\n",  idimension  ); 

/*  get  size  of  input  message  data  */ 
lwbgetparam(  in,  0,  LWBMSGSIZES,  insizes  ); 

if  (dx)  printf ("demod  insizes=%d,%d,%d\n",  insizes [0],  insizes [1], 

insizes [2]  ); 

count  =  insizes [0]; 

if  (  dx  )  printf (  "spewrit  count=%d\n",  count  ); 

/*  check  input  data  type  */ 
lwbgetparam(  in,  0,  LWBMSGTYPE,  Aitype  ); 
if  (  dx  )  printf (  "demod  itype=%d\n",  itype  ); 
if  (  itype  !=  5  ) 

{ 

lwberror(  ERROR,  "Writer  input  data  type  must  be  FLOAT"  ); 
return  (  0  ) ; 

} 

/*  open  output  files  */ 

if((outfilel  =  fopen(outf ilenaml , "a+"))  ==  NULL) 

{ 

lwberror(  ERROR,  "Couldn't  open  the  MAG  file"  ); 
return (  0  ) ; 

} 

if((outfile2  =  fopen(outf ilenam2,"a+"))  =  NULL) 

{ 

lwberror(  ERROR,  "Couldn't  open  the  FREQ  file"  ); 
return (  0  ) ; 

} 

if((outfile3  =  fopen(outfilenam3,"a+"))  =  NULL) 

{ 

lwberror(  ERROR,  "Couldn't  open  the  NOISE  file"  ); 
return (  0  ) ; 

} 

/*  release  the  input  */ 
lwbrelease(  in  ); 
re turn (  1  ) ; 

/*****»*****«**************  run  ************************************/ 

ufrun(  in  ) 
int  in; 

{ 

int  i,  bin; 

float  max  =0.0,  freq,  noise=0.0,  nave=0.0,  snr=0.0; 

if  (  ! lwbgetdatbuf  (  in,  Aidata  )  )  retum(  0  ); 

/*  process  data  */ 

for  (i=8179;  i  <  8208;  i++) 

{ 

if(idata[i]  >  max)  /*  get  max  amplitude  and  it's  freq  in  +/- 

5kts  band  about  57.0  Hz  */ 

{ 

max  =  idatafi] ; 
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} 


freq  =  (float) i  *  0.0069585; 
bin  =  i; 


} 

for  (i=6469;  i  <  9917;  i++) 

noise  =  noise  +  idata[i];  /*  get  noise  power  in  24  Hz  passband  */ 
nave=noise  /  3447; 

max  =  (max  /  1167.40902475);  /*  convert  from  LWB  to  reality  */ 

noise  =  (noise  /  1167.40902475); 
noise  =  noise  -  max; 
snr  =  max  /  noise; 

printf ("bin=%d  %fhz  m=%e  n=%4.4f  \n" .bin, freq, max, noise) ; 
printf ("  snr=%4.4f  ",snr); 

fprintf (outfilel ,"%4.4e  ”,max); 
fprintf (outf ile2,"%4.4f  ",freq); 
fprintf (outfile3,"%4.4e  ".noise) ; 
printf ("f ini shed  row  number  %d\n",m++); 
lwbrelease(  in  ); 
return (  1  ) ; 

} 

^* ****************************  STOP  *********************************/ 

ufstop(  in) 
int  in; 

{ 

if  (  ! in  )  retum(  0  ) ; 
lwbrelease(  in  ); 
printf ("\n"); 
fprintf (outfilel ,"\n") ; 
fprintf (outf ile2, "\n") ; 
fprintf (outfile3,"\n") ; 
fclose(outf ilel) ; 
fclose(outfile2) ; 
fclose(outf ile3) ; 

printf ("number  of  rows  =%d\n",m-l); 
return (  1  ) ; 

} 

!%% **************************  menu  **********************************/ 

ufmenu(  key  ) 
int  key; 

{ 

if  (  key  =  1  ) 

{ 

lwbgetstring(  "Enter  MAG  Filename",  1,  15,  outfilenaml  ); 
if  (  ! lwbupdatemenu(  1,  "%s",  outfilenaml)  )  retum(  0  ); 

> 

if  (  key  =  2  ) 

{ 

lwbgetstring(  "Enter  FREQ  Filename",  1,  15,  outfilenam2  ); 
if  (  ! lwbupdatemenu(  2,  "%s",  outfilenam2)  )  return!  0  ); 

} 

if  (  key  =  3  ) 
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ftthresh) ; 


lwbgetstring!  "Enter  NOISE  Filename",  1,  15,  outfilenam3  ); 
if  (  !  1  wbupdat emenu (  3,  "%s",  outfilenam3)  )  retum(  0  ); 

if  (  key  =  4  )  /*  not  currently  used  */ 

Iwbgetfloat!  "Enter  Output  Threshold",  1.0,  1000.0, 


if  (  ! lwbupdatemenu!  3,  "%f",  thresh  )  )  retum(  0  ); 


return!  1  ) ; 

/*******************************  drop  ********************************/ 

ufdropO 

{ 

return!  1  ) ; 
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